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A  New  Approach  for  the  Solution  of  Optimal  Control  Problems  on  Parallel  Machines 
Thesis  directed  by  Professor  Pieter  Frick 


This  thesis  develops  a  highly  parallel  solution  method  for  nonlinear  opti¬ 
mal  control  problems.  Balakrishnan’s  epsilon  method  is  used  in  conjunction  with  the 
Rayleigh-Ritz  method  to  convert  the  dynamic  optimization  of  the  optimal  control 
problem  into  a  static  optimization  problem.  Walsh  functions  and  orthogonal  poly¬ 
nomials  are  used  as  basis  functions  to  implement  the  Rayleigh-Ritz  method.  The 
resulting  static  optimization  problem  is  solved  using  matrix  operations  which  have 
well  defined  massively  parallel  solution  methods.  To  demonstrate  the  method,  a  vari¬ 
ety  of  nonlinear  optimal  control  problems  are  solved.  The  nonlinear  Raleigh  problem 
with  quadratic  cost  and  nonlinear  van  der  Pol  problem  with  quadratic  cost  and  ter¬ 
minal  constraints  on  the  states  are  solved  in  both  serial  and  parallel  on  an  eight 
processor  Intel  Hypercube.  The  solutions  using  both  Walsh  functions  and  Legendre 
polynomials  as  basis  functions  are  given.  In  addition  to  these  problems  which  are 
solved  in  parallel,  a  more  complex  nonlinear  minimum  time  optimal  control  problem 
and  nonlinear  optimal  control  problem  with  an  inequality  constraint  on  the  control 
are  solved.  Results  show  the  method  converges  quickly,  even  from  relatively  poor  ini¬ 
tial  guesses  for  the  nominal  trajectories. 'Parallel  results  show  the  method  also  offers 
significant  speedup  and  rapid  solution  times  even  on  the  small  hypercube  computer. 
Potential  parallelization  is  much  greater  than  eight  and  significantly  greater  speed-up 
might  be  realized  on  larger  parallel  machines. 
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CHAPTER  1 


INTRODUCTION 


The  theory  of  nonlinear  optimal  control  is  relatively  well  developed.  Un¬ 
fortunately,  due  to  the  high  computational  cost,  nonlinear  optimal  control  theory 
is  rarely  implemented  in  actual  control  hardware.  Instead,  nonlinear  systems  are 
linearized  and  linear  control  laws  are  used,  or  results  using  full  nonlinear  system 
equations  are  calculated  off  line  and  then  stored  in  memory  for  actual  implementa¬ 
tion.  Both  methods  result  in  suboptimal  forms  of  control.  In  recent  years,  parallel 
processors  have  become  well  developed.  The  fastest  computers  are  now  parallel  pro¬ 
cessing  computers.  As  parallel  processing  technology  develops,  it  should  offer  parallel 
processing  machines  with  enormous  computational  power  for  very  little  cost.  With 
inexpensive  parallel  processing  power  it  may  become  reasonable  to  consider  solving 
nonlinear  optimal  control  problems  in  real  time.  Real  time  solution  of  nonlinear  op¬ 
timal  control  problems  is  the  motivation  behind  developing  nonlinear  optimal  control 
methods  which  may  be  implemented  on  parallel  processing  machines. 

There  have  been  relatively  few  investigations  into  parallel  methods  of  solving 
optimal  control  problems.  This  is  probably  due  to  the  slow  development  of  parallel 
computers  and  the  lack  of  massively  parallel  methods  of  solving  the  system  dynamical 
equations.  The  following  is  a  review  of  parallel  methods  which  have  been  developed 
to  solve  optimal  control  problems 

The  first  research  into  parallel  processing  for  optimal  control  was  done  by 
Larsen  and  Tse  [1].  Their  method  was  to  give  a  parallel  implementation  of  dynamic 
programming  and  use  this  parallel  dynamic  programming  method  to  solve  the  non¬ 
linear  optimal  control  problem. 
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Travassos  and  Kaufman  [2]  chose  to  develop  a  parallel  implementation  of  the 
shooting  method  (variation  of  extremals).  They  used  parallel  initial  costate  searches 
along  with  a  parallel  state  equation  solver  to  formulate  a  parallel  method  for  solving 
nonlinear  two  point  boundary  value  problems.  They  also  proposed  the  type  of  parallel 
computer  architecture  needed  to  run  the  algorithm. 

Finally,  Meyer  and  Podrazik  [3]  looked  at  parallel  gradient  -  based  methods 
for  solving  the  discrete  optimal  control  problem.  They  used  parallel  gradient  search 
methods  along  with  a  parallel  linear  recurrence  solver  to  formulate  a  parallel  gradient 
method  for  discrete  systems. 

Two  major  issues  are  not  adequately  addressed  in  the  cited  literature;  in¬ 
terprocessor  communication  and  discretization  of  continuous  systems.  None  of  the 
above  researchers  actually  implemented  their  method  on  a  parallel  machine.  Because 
of  this,  the  very  important  issue  of  communications  between  processors  is  not  ad¬ 
equately  addressed.  Excessive  or  untimely  communication  between  processors  can 
severely  limit  the  speed-up  obtainable  on  a  parallel  machine.  Secondly,  all  the  above 
methods  require  that  the  system  be  discrete  or  the  system  be  simulated  in  some  dis¬ 
crete  manner  (state  equation  solving,  etc.).  For  some  very  nonlinear  problems  this  is 
a  significant  drawback  since  the  computational  load  to  obtain  a  discrete  system  from 
a  continuous  one  can  be  excessive.  This  dissertation  suggests  a  new  parallel  approach 
to  the  solution  of  the  optimal  control  problem  which  addresses  both  issues. 

The  fundamental  idea  of  the  new  method  is  simple.  Most  methods  for  solving 
two  point  boundary  value  problems  involve  solving  the  system  dynamical  equations 
many  times.  This  is  typically  the  most  computationally  intensive  part  of  these  meth¬ 
ods.  The  solution  of  ordinary  differential  equations  on  a  computer  is  inherently  serial. 
For  most  numerical  methods  of  solving  ordinary  differential  equations,  calculation  of 
the  next  point  is  dependent  on  current  points  or  past  points.  As  a  result,  solution  of 
the  dynamical  equations  introduces  a  type  of  serial  bottleneck  into  parallelizing  these 
types  of  methods.  In  order  to  obtain  a  solution  method  which  is  massively  parallel, 
the  dynamic  optimization  problem  is  converted  into  a  static  optimization  problem. 
This  conversion  is  obtained  by  use  of  Balakrishnan’s  epsilon  method  which  avoids 
solving  the  dynamical  equations.  The  new  static  problem  can  be  optimized  using 
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massively  parallel  methods.  Computational  results  are  very  promising.  The  method 
has  very  good  convergence  properties  and  has  highly  parallel  implementations. 

Chapter  2  develops  background  necessary  to  obtain  the  parallel  solutions 
for  optimal  control  problems.  Some  key  ideas  for  parallel  processing  for  numerical 
computation  are  presented,  Walsh  functions  are  reviewed,  and  an  introduction  to 
Balakrishnan’s  epsilon  method  is  given.  In  Chapter  3,  the  linear  optimal  control 
problem  with  quadratic  cost  is  formulated  as  the  epsilon  problem  and  is  solved  using 
the  Rayleigh  -  Ritz  method  with  Walsh  functions  as  a  basis.  This  conversion  gives 
a  problem  formulation  in  terms  of  matrix  operations  which  can  be  implemented  in 
parallel.  One  parallel  implementation  for  an  eight  node  Intel  Hypercube  is  given  and 
computational  results  for  a  linear  example  are  presented.  In  Chapter  4,  the  nonlinear 
problem  is  solved  by  quasilinearization.  Computational  results  for  a  nonlinear  optimal 
control  problem  with  quadratic  cost  are  given.  Also  in  chapter  4  the  method  is 
extended  to  include  the  case  of  terminal  constraints  on  the  state.  Chapter  5  reviews 
orthogonal  polynomials,  specifically  the  Legendre  polynomials,  and  gives  some  of 
their  properties.  The  Legendre  polynomials  are  then  used  as  basis  functions  to  solve 
the  two  nonlinear  examples  given  in  Chapters  4  and  5.  Chapter  6  shows  how  a 
variety  of  nonlinear  operations  can  be  performed  on  the  Walsh  coefficients  which 
approximate  a  function.  A  nonlinear  minimum  time  optimal  control  problem  is  then 
solved.  Chapter  7  uses  some  results  of  Chapter  6  to  solve  nonlinear  optimal  control 
problems  with  inequality  constraints.  A  numerical  example  is  given.  Chapter  8  draws 
final  conclusions  and  offers  some  additional  areas  of  research  which  may  be  of  interest. 
Appendix  A  reviews  some  systems  notation  used  throughout  the  thesis. 


CHAPTER  2 


BACKGROUND 


This  chapter  covers  the  background  necessary  to  develop  the  parallel  solution 
of  the  optimal  control  problem. 


Parallel  Processing  for  Numerical  Computation 


Much  research  has  been  done  in  the  area  of  parallel  processing  for  numerical 
computation.  In  this  section  an  approach  is  given  for  developing  parallel  implemen¬ 
tations  for  numerical  computation  problems.  Briefly,  the  idea  is  to  formulate  a  given 
numerical  problem  in  terms  of  smaller  operations  which  are  known  to  have  efficient, 
well  defined  parallel  implementations.  It  is  then  possible  to  combine  these  smaller 
operations  into  a  larger  parallel  scheme  which  may  retain  the  desirable  characteristics 
of  the  smaller  operations. 

There  are  three  properties  which  are  desirable  in  a  parallel  implementation, 
efficiency,  scalability,  and  flexibility.  Each  of  these  will  now  be  defined  and  addressed. 
Consider  first  some  performance  measures  for  parallel  implementations.  Let  T*(n ) 
be  the  time  to  solve  a  problem  serially,  where  n  is  some  measure  of  the  size  of  the 
problem.  For  example,  if  multiplying  two  nxn  matrices,  T*(n)  =  0(n3).  Let  Tp(n ) 
be  the  time  to  solve  a  problem  with  p  processors.  Then  the  speed-up  is  defined  as 

s'(n)  =  $j  <21> 


Ep(n) 


T*(n) 

PTp(n) 


SP(n) 

P 


And  the  efficiency  is 


(2.2) 
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so  Ep(n )  is  a  measure  of  the  time  a  processor  is  doing  work.  If  Ep(n )  is  bounded  away 
from  zero  as  n  and  p  increase,  then  as  n  increases,  p  can  be  increased  and  Sp(n)  will 
continue  to  increase.  For  example,  let 

Ep(n)  =  c 

where  c  is  some  constant  so  the  efficiency  is  bounded  away  from  zero.  Then 


Sp(n)  =  pc 

so  if  p  |,  then  Sp(n )  f.  If  c  J,  0,  efficiency  is  not  bounded  away  from  zero  and  speed-up 
may  not  increase.  If  efficiency  is  bounded  away  from  zero  as  n  is  increased  without 
bound,  Sp(n )  will  continue  to  increase  indicating  that  the  number  of  processors  which 
can  be  used  to  solve  the  problem  can  be  increased  without  bound.  If  efficiency  goes  to 
zero,  Sp(ra)  may  not  increase  indicating  that  adding  more  processors  may  not  speed 
up  solution  and  therefore  might  be  a  wasteful  endeavor.  A  problem  which  can  be 
solved  in  parallel  and  has  efficiency  bounded  away  from  zero  will  be  called  scalable. 

One  very  important  characteristic  of  parallel  computers  which  effects  scala¬ 
bility  is  communication  time.  Communication  time  is  the  time  to  send  one  message, 
typically  a  real  number,  from  one  processor  to  another.  It  is  very  important  to  con¬ 
sider  communication  time  when  the  solution  time  of  a  particular  problem  is  of  interest. 
A  parallel  solution  which  demands  the  passing  of  many  messages  can  quickly  become 
communication  bound.  A  problem  which  is  communication  bound  spends  most  of 
the  time  communicating  and  little  time  computing.  Extra  time  spent  in  excessive 
communications  can  quickly  render  a  problem  unscalable.  Therefore,  a  problem  will 
only  be  considered  scalable  if  it  is  scalable  including  communication  time. 

Finally,  consider  the  parallel  architecture  which  is  necessary  to  run  a  partic¬ 
ular  numerical  problem.  Very  often,  solution  of  a  particular  problem  is  bound  strictly 
to  a  particular  type  of  parallel  architecture.  A  problem  which  can  be  solved  on  a  va¬ 
riety  of  parallel  architectures  will  be  said  to  have  a  flexible  parallel  implementation. 


Problem 

Time 

Topology 

Inner  Product 

Matrix- Vector  Multiplication 
Matrix-Matrix  Multiplication 
Q  R  Factorization 
Backsubstitution 

0(log  n) 
0(n) 

0(n) 

0(n) 

0(n) 

Hypercube  w/  p=n  processors 
Linear  Array  w/  p=n  processors 
Mesh  w /  p  =  nJ  processors 

Mesh  w/  p  =  n2  processors 

Linear  Array  w/  p=n  processors 

Table  2.1:  Times  for  Parallel  Matrix  Calculations 


It  is  obviously  most  desirable  for  a  problem  to  have  a  flexible  parallel  implementation 
since  then  the  problem  may  be  solved  on  a  variety  of  parallel  computers  or  parallel 
processing  arrays.  The  goal  then  is  to  strive  for  a  parallel  implementation  which  is 
efficient,  flexible,  and  scalable. 

Many  simple  matrix  operations  have  the  desirable  properties  of  efficiency, 
flexibility,  and  scalability.  For  example,  matrix  multiplication  and  the  solution  of  sys¬ 
tems  of  linear  equations  have  a  variety  of  implementations  both  on  parallel  computers 
and  array  processors  such  as  systolic  arrays  or  wavefront  arrays  and  have  been  shown 
to  be  scalable  [4,  5].  Table  2.1  (adapted  from  [4])  gives  some  examples  of  other  sim¬ 
ple  matrix  operations  which  are  known  to  have  efficient,  flexible,  and  scalable  parallel 
implementations. 

The  idea  then  is  to  break  down  solution  of  the  optimal  control  problem  into 
smaller  operations  which  have  well  defined  parallel  implementations  such  as  simple 
matrix  operations.  In  this  way,  it  is  possible  to  form  parallel  implementations  of 
optimal  control  problems  which  may  retain  the  desirable  parallel  characteristics  of 
the  simple  matrix  operations. 

Walsh  Functions  for  Optimal  Control 

Walsh  functions  have  been  used  for  systems  analysis,  systems  identification, 
and  optimal  control  [6,  7,  8].  The  following  is  a  review  of  some  of  their  useful  prop¬ 
erties  for  use  in  the  parallel  solution  of  optimal  control  problems. 

The  Walsh  functions  are  a  set  of  square  wave  functions  which  are  orthogonal. 
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Figure  2.1:  Rademacher  Functions 

The  Walsh  functions  axe  generated  from  a  more  elementary  set  of  square  waves  known 
as  the  Rademacher  functions  [9] .  Rademacher  functions  are  a  set  of  square  waves  of 
unit  height  which  have  periods  of  1,  , 21-fc.  The  first  four  Rademacher  func¬ 

tions  are  shown  in  Figure  2.1.  Rademacher  functions  have  odd  symmetry  about  t  =  0 
and  t  =  |.  Asa  result,  the  functions  are  incomplete  since  it  is  not  possible  to  approx¬ 
imate  a  function  which  has  even  symmetry  about  t  =  0  and  t  =  |  using  Rademacher 
functions.  Walsh  combined  the  Rademacher  functions  to  form  a  complete,  orthogonal 
set  of  rectangular  waves  [10]. 

The  Walsh  functions  can  be  generated  as  the  product  of  n  Rademacher 
functions.  Let  a„  •  •  •  a2ax  be  an  n  digit  binary  number,  so  the  a,’s  are  either  0  or  1 
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Figure  2.2:  Walsh  Functions 
and  an  is  1.  The  Walsh  functions  are  obtained  by 

&„...«,«»( 0  =  [r„ (<)]“”  •  •  •  [ra (OPM*)]*1  (2.3) 

EXAMPLE: 

4>o(t)  =  To(t)  </>10o(f)  =  **3(t) 

&(0  =  ri(0  </>ioi(0  =  i"3(0ri(0 

^io(0  =  r2(t)  <t>uo(i)  =  r3(f)r2(f) 

<£n(0  =  r2(t)ri(t)  ^ni(t)  =  r3(i)r2(t)ri(t) 


Figure  2.2  shows  the  first  four  Walsh  functions  <f>0(t )  to  $3(f).  Just  as  a  function  /(f) 
may  be  expanded  into  a  Fourier  series,  a  function  /(f)  which  is  integrable  in  (0,1] 


9 


may  be  expanded  into  a  Walsh  series. 

/(f)  =  co<f>o(t)  4-  ci<f>x{t)  +  •  •  •  +  CN-\4>N-\(t)  +  •  •  •  (2-4) 

Walsh  functions  have  very  useful  integration  properties.  Integration  of  a 
vector  of  Walsh  functions  can  be  performed  by  a  matrix  multiplication,  specifically 


f  y(r)dT  =  PV(t) 
Jo 


where 


tf(r)  = 


Mt ) 


\  <j>N-i(r)  ) 

P  is  called  the  integration  matrix  and  has  the  following  general  structure: 

_  i  pN  _  Pff/2  -(dv)7Af/2  N  =  2n 

~  2  *  l(^)/jyr/2  0  J  n  =  1,2,... 

EXAMPLE:  Let  N  =  A  and  /(f)  =  1.  Find 


/  f(r)di 

Jo 


Then 


I  -i  -i  o 

P=  i  0  0  -I 

i  o  o  0 


0^00 


Let  the  Walsh  approximation  for  /(f)  be 


where 


m  =  Fv(f) 


F=  [  1  0  0  0 


[*  FV(T)dT  =  F  ['  ¥(r)dT  =  FP*(t) 
Jo  Jo 


then 


10 


where 


"=[§  -}  -J  o] 

and  it  can  be  shown  that  F P^{t)  is  the  four  term  Walsh  series  approximation  for 

g(t)  =  t- 


Walsh  functions  have  interesting  multiplication  properties.  The  product  of 
two  Walsh  functions  is  found  by  the  mod  2  addition  of  the  binary  form  of  the  Walsh 
function  subscripts.  Thus  the  Walsh  functions  form  a  closed  set  under  multiplication. 
The  multiplication  of  two  Walsh  functions  never  results  in  a  Walsh  function  whose 
subscript  is  greater  than  the  functions  being  multiplied. 


EXAMPLE: 


<f>io(t)  x  <f>n(i)  —  <t>o\ (t) 

Finally,  the  Walsh  functions  form  an  orthogonal  system  on  0  <  t  <  1 ,  or 

0  i  ±  j 


f1  0  1*7 

/  j(i)A  =  \ 

Jo  I  1  i-j 


(2.8) 


The  Epsilon  Method 

In  order  to  introduce  the  general  mathematical  structure  for  using  the  ep¬ 
silon  method  consider  the  following  general  optimal  control  problem.  It  is  desired  to 
control  some  system  described  by: 

x(t)  =  f(x{t),u(t),t)  x(t0)  =  x  o  (2.9) 

where  x{t )  6  Rn,u{t )  €  R"1  and  u(t)  represents  the  control.  A  cost  functional  V 
describes  some  performance  measure  to  be  met.  In  general, 

V  =  [*'  G{x(t),u(t),t)dt  (2.10) 

Jo 

and 

v>(  x(t,),<,)  =  0  (2.11) 


is  a  terminal  constraint  on  the  state  which  must  be  met  exactly.  The  problem  is  to 
find  u*(t)  which  minimizes  the  cost  function  while  meeting  the  dynamic  constraint  of 

(2.9). 
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A  commonly  used  approach  to  solve  the  continuous  optimal  control  problem 
is  to  employ  the  method  of  Lagrange  multipliers  [12]  to  adjoin  (2.9)  and  (2.11)  to  the 
cost  functional 

V'  =  if)+  (  1  [G(x(t),  u(t),  t)  +  XT(t){f(x(i),  u (t),  t)  -  x(t)}]dt  (2  12) 

Jo 

Since  (2.9)  must  hold  at  each  t  £  [0,  tf]  the  multiplier  X{t)  £  BP  is  a  function  of  time 
and  is  known  as  the  costate.  Equation  (2.11)  holds  only  at  one  time  and  so  u  is  a 
constant  multiplier. 

Now  define  the  Hamiltonian  function  as 


H(x(t),u(t),\(t),t)  =  G(x(t),u(t),t)  +  XT(t)[f{x(t),u(t),t)]  (2.13) 


Then 

V'  =  vTi>(x(tf),tf)  -f  f  /[if(x(f)>u(t),A(t),t)-  XT(t)x(t)]dt  (2.14) 

J  0 

Now  minimizing  V  subject  to  the  constraint  of  equation  (2.9)  is  equivalent  to  the 
unconstrained  optimization  problem  of  minimizing  V' .  It  can  be  shown  that  the 
following  conditions  must  be  met  to  obtain  a  minimum  of  V' . 


m  =  aj/(l(‘)j“(])),A(f),‘)  =  /(*«.<* <o.o  t  >  o 

-x(t)  =  dH{xWg®’Kt)’i)  t<t, 


dH(x(t),u(t),X(t),t) 

du(t) 

x(0)  =  Xo 


=  0 


(2.15) 

(2.16) 
(2.17) 


O/'Ji'  -  A(0)T)  It;  dl(*/)  +  (V’J*'  -  H(x(i),u(t),X(t),t))  I tf  dtf  =  0  (2.18) 


In  (2.18)  the  subscripts  denote  partials  and  the  bar  means  the  quantity  is  evaluated 
at  tf. 

The  continuous  optimal  control  problem  has  been  converted  into  a  two  point 
boundary  value  problem  (TPBVP).  The  state  and  costate  form  a  system  of  differential 
equations  in  which  the  initial  conditions  on  the  state  and  the  final  conditions  on 
the  costate  are  known.  The  structure  of  the  TPBVP  suggests  solution  methods 
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which  involve  numerically  solving  the  state  and  costate  equations  many  times.  An 
example  of  this  is  given  in  Chapter  6.  Since  numerical  methods  for  solving  ordinary 
differential  equations  are  very  difficult  to  parallelize,  this  introduces  a  serial  bottleneck 
into  parallelization  of  the  problem.  A  method  which  avoids  solving  the  state  equations 
was  formulated  by  Balakrishnan  [13]. 

Balakrishnan’s  epsilon  method  treats  the  dynamic  equations  as  a  penalty 
equality  constraint.  More  specifically,  in  the  optimal  control  problem  it  is  desired  to 
minimize  some  cost  function 

V(x(t),u(t))  =  f  '  G(x(t),u(t),tyt  (2.19) 

Jo 

while  meeting  the  constraints  imposed  by  some  dynamical  system 

x{t )  =  f(x(t),u(t),i)dt,  x(0)  =  x,,  (2.20) 

where  t  €  [0  <  oo.  By  treating  the  dynamical  equations  as  a  penalty  constraint, 

the  problem  can  be  reformulated  into  a  nondynamic  form: 

J(x(t),u(t)\e)  =  I!  ||2  dt  +  V(x(t),u(t))  (2.21) 

where 

e(t;e)  =  x(t)  -  f(x(t),u(t),t)  (2-22) 

The  composite  cost  functional  (2.21)  is  minimized  with  respect  to  x(t)  and  u(t)  for 
a  given  e  >  0.  If  necessary,  the  process  is  performed  with  a  sequence  {e^}  which 
decreases  monotonically.  The  process  is  stopped  when  the  error  in  the  system  dy¬ 
namics  is  small.  This  formulation  has  been  called  the  differential  epsilon  method. 
Balakrishnan  made  a  thorough  study  of  the  differential  epsilon  method  and  gave  the 
conditions  where  the  solution  to  the  epsilon  problem  approaches  the  optimal  solution 
to  the  optimal  control  problem  as  c  — >  0  [13]. 

An  integral  formulation  of  the  epsilon  method  was  developed  by  Frick  [14] 
By  adjoining  an  error  function  described  by 

e(t-,e)  =  x(t)  -  x.  -  f  /(x(r ),  u(r ),  t )dr 
Jo 


(2.23) 
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to  the  cost  function,  the  composite  cost  functional  is  obtained: 

As, *(<),»(!))  =  £'  II  ll!  dt  +  (2.24) 

This  composite  cost  function  is  minimized  in  the  same  fashion  as  the  differential 
epsilon  method.  Convergence  for  linear  optimal  control  problems  with  quadratic 
performance  indices  was  shown  by  Frick  [14].  These  results  are  summarized  next. 
Convergence  Result: 

Consider  the  sequence  of  scalars  {£j}  ].  0  monotonically  decreasing.  For  the 
corresponding  sequence  of  minimizing  solutions  to  the  integral  e  problem  denoted  by 
{x0(£j),  u0(£j)}  and  the  associated  {e0(£j)}  >  we  have: 

u°(ej)  u*  while  xo{Cj)  x*  and  °(  — >  A*  (2.25) 

cj 

e0{Cj)  |  0  while  J(evx0(ej),u0(e3))  T  V(x*,u*)  (2.26) 

as  £  J,  0  .  Here  x*,u*  and  V(x*,u*)  are  the  optimal  control,  state,  and  cost  for  the 
case  of  linear  dynamical  equations  and  quadratic  cost. 

For  practical  application,  the  Rayleigh-Ritz  method  has  been  used  to  min¬ 
imize  the  composite  cost  functional  of  the  differential  epsilon  method  [18].  Gradient 
techniques  have  been  used  to  minimize  the  composite  cost  functional  of  the  integral 
epsilon  method  [14].  In  the  next  chapter,  the  integral  epsilon  problem  for  the  case  of 
linear  time- varying  dynamical  equations  and  quadratic  cost  function  is  solved  using 
the  Rayleigh-Ritz  method  with  Walsh  functions  as  basis  functions.  As  a  result  of 
the  integration  properties  of  the  Walsh  functions,  the  solution  of  the  integral  epsilon 
problem  using  the  Rayleigh-Ritz  method  is  an  approach  which  allows  for  a  high  degree 
of  parallelism  in  the  solution. 


CHAPTER  3 


PARALLEL  SOLUTION  OF  LINEAR  OPTIMAL  CONTROL 

PROBLEMS 


In  this  chapter,  the  epsilon  problem  for  the  case  of  linear  time-varying  dy¬ 
namics  and  quadratic  cost  is  solved  by  the  Ritz  method  using  a  Walsh  function  basis 
set.  The  resulting  static  optimization  problem  can  be  solved  using  matrix  operations. 
A  parallel  implementation  for  a  linear  array  embedded  in  a  hypercube  computer  is 
then  developed. 


The  Rayleigh-Ritz  Method  using  a  Walsh  Basis 

To  minimize  a  cost  functional  using  the  Rayleigh-Ritz  method,  the  time- 
varying  functions  in  the  cost  functional  are  approximated  by  a  series  of  basis  functions. 
Approximation  by  a  Walsh  series  will  be  used.  Once  the  time-varying  functions 
in  the  cost  functional  have  been  replaced  by  the  parameterized  functions,  the  new 
cost  functional  can  be  minimized  with  respect  to  the  constants  of  the  approximated 
functions.  In  this  manner  minimization  of  the  cost  functional  is  replaced  by  a  static 
optimization  problem  which  is  very  amenable  to  parallel  implementation. 

It  is  desired  to  control  some  dynamical  system  which  is  described  by  a  set 
of  linear  time  varying  differential  equations. 

x(t)  =  A{t)x{t)  +  B(t)u(t)  +  C(t),  x(0)  =  x.,  (3.1) 

with  x(t)  €  G  fl™.  A  performance  index  or  “cost  function”  is  selected  to 

force  the  dynamical  system  to  follow  some  desired  performance  specification.  For 
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example, 

V(x(t),  u(i))  =  z(0.  Q*(0  >  <  “(<).  Ruit)  >}dt  (3.2) 


Here  <  x(t),  Qx(t)  >  denotes  the  inner  product.  Q  is  positive  semi-definite  and  R 
is  positive  definite.  This  is  a  quadratic  cost  function.  It  penalizes  deviations  from 
zero  in  the  state  and  the  control,  especially  if  Q  and  R  have  large  magnitude.  The 
linear  time-varying  optimal  control  problem  is  to  find  the  optimal  control  u*(f )  which 
drives  the  states  of  the  dynamical  system  along  a  trajectory  x*(t)  such  that  the  cost 
functional  is  minimized  and  the  dynamical  equations  are  met.  This  linear  time- 
varying  optimal  control  problem  can  be  solved  using  the  integral  epsilon  method  and 
convergence  of  the  method  is  assured.  Consider  the  integral  form  of  (3.1) 


x(t)  =  x,  +  /  .A(t)x(tWt  +  f  B(r)u(r)dr  - f  /  C(t)cLt  (3.3) 

Jo  Jo  Jo 


These  dynamical  equations  can  be  viewed  as  a  dynamical  constraint  and  can  be 
adjoined  to  the  cost  functional  by  forming  an  error  function 


;(f;  e)  =  x(t)  —  x,  —  f  v4(t)x(t)<£t  —  f  B(t)u(t)cIt  —  /  C(r)dr  (3.4) 

Jo  Jo  Jo 


Then  the  composite  cost  functional  can  be  formed 

J(c,x(t),u(t))  =  Jo'{^  II  e(i;£)  f  +i  <  x(t),Qx(t)  >  <  "(*).*<<)  >}* 

(3,5) 

The  first  step  is  to  expand  the  time  varying  functions  in  the  composite  cost 
function  into  a  Walsh  series.  Since  for  computational  application  the  series  must 
be  truncated,  the  Walsh  series  represents  an  approximation  to  the  original  function. 
Once  the  problem  is  in  approximation  form,  it  is  solved  by  finding  the  coefficients 
to  the  state  and  control  which  minimize  the  composite  cost  function.  The  Walsh 
functions  form  an  orthogonal  base  on  [0,1)  so  assume  (without  loss  of  generality)  that 
tf  —  1.  Now  approximate  x(t),u(t),  A(t),  and  B(t)  with  N  Walsh  functions  where 
N  =  2k,k  >  0.  For  example: 


z(0  —  co<t>o(t)  +  Ci<f>i(t)  +  •  •  •  +  CS-l4>N-l(t) 
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The  following  approximations  can  be  made 


x{t)  «  X9(t) 
u(t)  »  U*(t) 
A{t)  «  a^(i) 
S(t) 


(3-6) 


where  X,  17,  a,  and  /3  axe  matrices  of  coefficients  and  ^(t)  =  [<£0(f),  . . .  <j>N-\{t)}T  ■ 

Also,  the  error  term  in  the  composite  cost  functional  can  be  approximated: 


e(£;  e )  «  Etyft) 


(3.7) 


Equation  (3.5)  can  be  written  in  approximation  form: 

J=  /1{^-^T(£)f;rJE;^(£)  +  ^T(i)A:T<?X^(£)  +  ^T(£)t/Ti217^(£)}d£  (3.8) 

J  o  2e  2  2 


The  matrix  product  terms  have  the  form  ^T(t)M ^(T)  which  can  be  written  in  sum¬ 
mation  notation 

N-l  N-l 

(3.9) 


yT{t)MV(t)  =  &(0  E  rriij<f>j{t) 

»= 0  j=0 


and  therefore 


['  tfT(£)M¥(£)df  =  ZE  mv  f 1  Ut)<l>i{t)dt  (3.10) 

*'0  j=0  j=0  "'O 

As  noted  in  Chapter  2,  for  the  Walsh  functions  the  above  integral  is  the  identity 
matrix  and  so 

.1  N-l 

/  'JfT(£)A/^r(£)d£  =  ^2  mu  =  trace(M)  (3.11) 

Jo  ._n 

Using  (3.11),  (3.8)  can  be  written 


»=o 


J  =  trace{EETE  +  \xTQX  +  \uT  RU) 


(3.12) 


Now  it  is  possible  to  take  advantage  of  some  powerful  systems  notation  standardized 
by  Brewer  [15]  (See  Appendix  A).  Using  vec  notation  (3.12)  can  be  written 


J  =  —vec(E)Tvec(E)  +  ^-vec(X)T  vec{QX)  +  \-vec{U)Tvec{RU) 

X/C  Z  Z 


(3.13) 
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or 

J  =  vec{E)Tvec(E)  +  ]rvec(X)T(IN®Q)vec{X)+}-vec(U)T(IN®R)vec(U )  (3.14) 

It  is  desired  to  find  the  Walsh  coefficients  of  the  state  and  control  vec*(X)  and  vec*(U) 
which  minimize  the  above  cost  function.  Since  vec(E)  is  a  function  of  vec(X)  and 
vec(U),  it  must  be  found  next. 

To  approximate  e(t;  e)  in  terms  of  Walsh  functions,  the  time- varying  in¬ 
tegrals  in  (3.4)  must  be  approximated.  This  can  be  done  using  the  multiplication 
and  integral  properties  of  Walsh  functions.  Consider  first  the  scalar  case  where  n=l. 
Let  A{t)  «  a’I'(f)  or  tyTaT,  a  is  a  vector  of  coefficients,  a  =  [a0,  <*i,  •  •  •  ow-i]  and 
x(f)  xs  A'^'(t)  Then 

A(t)x(t)  «  Xtf(f)tfT(f)aT  (3.15) 

Using  the  rule  given  earlier  for  multiplying  Walsh  functions,  the  matrix 
can  be  found.  For  example  if  iV  =  4 

M *)  M*)  Mi) 

*w(o*r.)(o  =  m  Mt)  m  m 

1  ’  ( ’  hw  hw  hw  hw 

_  Mi)  Mi)  Mi)  MO  . 

Chen  [6]  gives  the  general  form  of  the  product  of  these  two  matrices  as 


*<?)+¥« 

*(f) +?«<•?*,+*(<)  *(f) 


(3.16) 


where  the  subscript  (|)+  y  denotes  that  the  subscript  of  each  element  is  increased 
by  y.  Notice  that  the  elements  of  this  matrix  have  a  single  subscript.  For  matrices 
of  this  type,  the  following  is  true 


¥(i)«T(t)a  =  A  a*(t) 


(3.17) 


Here  a  is  some  vector  of  coefficients  and  Aa  is  a  matrix  of  the  elements  of  a  which 
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have  the  same  subscripts  as  ^r(t)^,T(t).  For  example  if  N  =  4: 


4>o{t)  fa  (t)  fa{t)  fa(t) 

Qo 

a0fa(t)  4  Qifa(t)  4  Q-ifaii)  4  q3^3(0 

fa  (t)  <t>o(t)  fa (0  fa{t) 

Otl 

oi\fa{t)  -f  aofa(t)  +  a $fa{t)  -(-  a.2fa{t) 

fa  (i)  fa(t)  fa(t)  fa  (t) 

a  2 

Ot  2fa(t)  4  <*3^1  (0  +  OLofa(t)  +  <>1^3(0 

fa(t)  fa (t)  fa(t)  fa (*)  _ 

03 

Q!3^o(0  4  oc2<f>i(t)  4  aifa{t)  4-  &o<f>3 (t) 

and 


<*ofa(t)  4  ai^i(0  4  O-ifa (t)  4  <23^3 (0 

ao  a.\  a2  0:3 

fa(t) 

&ifa{t)  4  cxo<f>i(t )  -f  a3<f>2[t)  4  ct2fa(t) 

Qi  ctQ  CL3  a  2 

fa{t) 

ot2<f>o(i)  4  ot3<f>i(t)  +  ao<f>2(t )  4  ct\fa(t) 

Ct2  CL  3  a0  CL\ 

fait) 

03^0(0  4  a2fa(t)  +  ai<^2(i)  4  aofa(t) 

Ct3  a2  0(1  O(0 

_  &(0 

Now  (3.15)  can  be  rewritten 

A(t)x(t)*z  XAaV(t)  (3.18) 


then  for  the  scalar  case 

f  A{t)x{t)<1t  tsXAaP*(t) 
Jo 

and 

vec(XAaP)  =  [(AaP)T]uec(X) 
Now  consider  the  case  where  n  >  1,  then 


an(<)  a12(t)  •••  aln(t) 

Xl(t) 

A(i)  = 

a2j(t) 

;  *(0  = 

x2(t) 

ayni(t)  ■  ■  •  ann(f) 

Xn(0  . 

Using  (3.19),  the  time-varying  integral  for  n  >  1  can  be  written  as 

X1AanP^(t)-{-  •••  +XnA  ainP*(t) 

*,Aa„  ?#(<)+  +X„K„P*(t) 


(3,19) 


(3.20) 


(3,21) 


(3.22) 
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where  an  =  [aj1,aj1,  •  •  •  aff1].  Here  the  superscript  denotes  the  Walsh  function 
number  so  au(i)xi(f)  =  a'n,Ir(t),^T^i-  Now  $(<)  is  integrated  out  of  (3.22)  as  shown 
in  (3.11)  when  the  cost  function  was  approximated  with  Walsh  functions.  Equation 
(3.22)  can  now  be  written, 

uec(-Xi) 

[(/„  ®  Pt)A1)  V‘C(-X2)  (3.23) 

vec(Xn) 


where 


and  this  can,  in  turn,  be  written  in  vec  notation: 


(3.24) 


vec(X/i (,““>P)  =  \(PT  8  7„)A<1M‘”)]vec(.X’)  (3.25) 


where  (if  N=4) 

'  AS  AJ  AS  AS  ’ 

A  _  AJ  AS  AS  AS 
“  AS  AS  AS  AS 
A l  A l  AI  A° 

Here  A°  denotes  a  matrix  of  the  zeroth  Walsh  coefficients.  The  order  of  subscripts  of 
Aq  is  the  same  as  the  scalar  case.  Now  vec(E)  can  be  written: 

vec(E)  =  {vec(X)  -  vec(G,)  -  \(PT  0  Jn)AQ]vec(A:) 

~[(PT  0  In)A0]vec(U)  -  ( PT  0  In)vec(C)} 

To  minimize  J  the  gradient  is  calculated  simultaneously  with  respect  to  both  vec(X) 
and  vec(U )  and  set  equal  to  zero. 


VuecWy  =  i(/wn  -  ( PT  <0  In)Aa)T  {(/Nn  ~  (PT  0  In)Aa)vec(X) 

~([PT  <0  In]A0)vec(U)  -  uec(G,)|  +  (Is  0  Q)vec(X)  =  Q 


(3.27) 
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Vvec([/)j  =  -\{PT  ®  7n] AP)T  {(Isn  -  ( PT  ®  In)Aa)vec(X) 

-([PT  ®  7n]A^ )vec(U)  -  vec(G,)}  +  (7^  ®  R)vec(U)  =  0 
The  resulting  system  of  equations  can  be  written  in  very  simple  form: 

'  tfu  Tfjj 

K21  K  22 

where 

*11  =  WNn  ~  (PT  ®  7n)Aa)T(/*n  ~  (PT  ®  7n)Aa)  +  (7*  ®  Q) 

*12  =  -i(7^„  -  (PT  ®  7n)Aa)T([PT  ®  JJA*) 

*21  =  -\([PT  ®  7„]A^)T(7yn  -  (PT  0  7n)Aa) 

(3.30) 

*22  =  i([PT  ®  7n]A^)T([PT  ®  7„]A^)  +  (Is  ®  P) 

7?r  =  i(7Nn  -  (PT  ®  7n)Aa)T(vec(G.)  +  [PT  ®  7n]vec(C)) 

D2  =  -i([PT  ®  7„]A^)T(uec(G,)  +  [ PT  ®  7n]vec(C)) 

Here  the  subscripts  on  the  identity  matrices  indicates  their  dimensions.  K  can  be 
written  more  concisely  as  follows,  if 


vec*(X) 

vec*(U) 


-  S 


(3.29) 


and 


Then 


7jvn  —  (PT  ®  7n)Aa  -(PT  ®  7n)A^ 
0  0 

L=  \  1 n®Q  0 

0  7*  ®  R 

K  =  -MtM  -+  7/ 

t 


(3.31) 

(3.32) 


(3.33) 


The  linear  time  varying  optimal  control  problem  has  been  reduced  to  Kronecker  prod¬ 
ucts,  matrix  multiplications,  and  the  solution  to  a  system  of  linear  algebraic  equations. 


21 


As  mentioned  in  Chapter  2,  these  types  of  operations  have  desirable  properties  for 
parallel  implementations.  To  conclude  this  section  the  method  is  illustrated  by  a 
simple  numerical  example. 

EXAMPLE: 

Solve  the  following  optimal  control  problem  using  the  epsilon  method  and 
the  first  4  Walsh  functions  as  basis  functions. 

x  =  x  -f  u  x(0)  =  1.0 

V(x,u)  =  f  x2  +  u3dt 
Jo 

So  N  =  4  and  the  following  is  true 

§  -i  -I  0  ' 

}  0  0  -i 

I  0  0  0 

0  I  0  0 

a=[l  0  0  o]  /3=[l  0  0  o]  c  =  [  1  0  0  o] 

Q  =  1  R  =  1  Aq  =  U  A  P  =  I4 

1 
0 
0 
0 

Thus  the  following  is  true, 

Ian  -  ( PT  ®  /n)Aa  =h-PT 

and 

-(PT®In)Ap  =  -PT 

Finally, 

U-PT  -PT' 

0  0  ~ 


22 


0.5000 

-0.2500 

-0.1250 

0.0000 

0.2500 

1.0000 

0.0000 

-0.1250 

0.1250 

0.0000 

1.0000 

0.0000 

0.0000 

0.1250 

0.0000 

1.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-0.5000 

-0.2500 

-0.1250 

0.0000 

0.2500 

0.0000 

0.0000 

-0.1250 

0.1250 

0.0000 

0.0000 

0.0000 

0.0000 

0.1250 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

h  0 
0  u 


L  = 

Now  choose  epsilon  small  to  ensure  small  dynamical  equation  error  and  calculate  K , 


K 

1  T 

=  -Mt 

s; 

4- 

ti 

II 

1 

£ 

.001 
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Figure  3.1:  Numerical  Example  -  Four  Walsh  Functions 


Solving  the  system  of 


equations  which  have  been  formulated  gives 


vec(X)  = 


0.9638 

-0.0549 

-0.0283 

0.0584 


-0.7387 

-0.4127 

-0.2128 

-0.0447 


The  resulting  Walsh  approximations  for  the  state  and  control  are  plotted  against  the 
continuous  solutions  in  Figure  3.1. 

In  the  next  section,  a  parallel  implementation  is  given  for  the  solution  of 
optimal  control  problems  on  an  Intel  Hypercube  with  eight  processors. 
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Parallel  Implementation 

In  this  section,  a  specific  parallel  implementation  for  solution  of  optimal 
control  problems  on  an  eight  processor  Intel  Hypercube  is  developed.  This  is  done 
by  first  considering  parallel  solution  of  each  of  the  necessary  matrix  operations  on  a 
linear  array  of  processors.  Once  this  is  done,  it  is  a  straightforward  task  to  develop 
a  specific  Hypercube  implementation  since  a  set  of  processors  connected  in  a  linear 
array  is  a  subset  of  a  set  of  processors  connected  in  a  hypercube  configuration. 

The  matrix  operations  of  interest  are  Kronecker  products,  matrix  multipli¬ 
cation,  and  the  solution  to  a  set  of  linear  algebraic  equations.  To  solve  the  set  of 
equations,  QR  triangularization  and  backsubstitution  are  used  because  of  the  stable 
numerical  properties  of  QR  triangularization.  The  serial  time  to  perform  a  matrix 
multiplication  of  two  n  x  n  matrices  and  the  serial  time  to  perform  a  QR  triangular¬ 
ization  on  an  n  x  n  matrix  are  the  same,  T*(n3).  Serial  time  for  a  backsubstitution 
is  T*(n2).  On  a  linear  array  of  processors,  matrix  multiplication  and  QR  triangular¬ 
ization  have  execution  times  of  T*(n2).  Since  the  parallel  QR  triangularization  and 
matrix  multiplication  have  solution  times  on  the  same  order  as  serial  backsubstitution, 
it  is  not  necessary  to  reduce  execution  time  of  the  backsubstitution  by  performing  it 
in  parallel.  It  may  be  performed  serially  and  not  significantly  effect  overall  execution 
time.  The  Kronecker  product  for  two  n  x  n  matrices  is  an  T*(n4)  time  operation.  The 
matrices  used  in  the  Kronecker  products  for  solution  of  the  optimal  control  problem 
are  typically  sparse.  This  makes  performing  the  full  Kronecker  products  very  ineffi¬ 
cient.  Consequently,  for  an  efficient  implementation  on  an  eight  node  Hypercube ,  the 
Kronecker  products  are  done  in  serial  since  they  are  not  computationally  intensive 
enough  to  justify  parallel  execution.  The  main  computational  burden  of  the  problem 
lies  in  the  matrix  multiplication  and  QR  triangularization.  Parallel  execution  of  these 
operations  is  considered  next. 
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Parallel  Matrix  Multiplication 


It  is  desired  to  do  some  matrix  multiplication  AB  =  C .  Note  that  B  can  be 
partitioned  by  columns 


B  —  [Bi  B2  •  ■  ■  Bn 


(3.34) 


Bi,  B2, . . .  Bn  are  the  columns  of  B  or  groups  of  columns.  The  columns  of  C  can  be 
found  as  the  product  of  A  and  columns  of  B. 


C  =  [ AB j  AB2  . . .  ABn ]  (3.35) 

These  matrix  multiplications  are  independent  and  can  be  done  in  parallel.  The  matrix 
A  and  appropriate  column  or  group  of  columns  are  passed  to  a  processor  which 
calculates  the  column  or  group  of  columns  of  C.  If  B  has  n  columns,  then  up  to  n 
processors  can  be  used. 


Parallel  Solution  to  Linear  System  of  Equations 


Typical  noniterative  methods  for  solving  systems  of  equations  Ax  —  b  consist 
of  triangularizing  the  augmented  matrix  [^4|fe]  and  using  backsubstitution  to  find  x. 
Bojanczyk  gave  a  parallel  implementation  for  QR  triangularization  on  a  2  dimensional 
mesh  of  processors  [16].  This  implementation  can  be  modified  for  a  linear  array  of 
processors. 

Assume 


an 

a2i 

an 

a»+i,l 


a12 

<*22 

Oi2 

o,+:,2 


au 

° 2} 

an 

o.+lj 


fllp 

a2p 


and  is  nonsingular.  An  orthogonal  matrix  Q  will  rotate  A  to  a  triangular  matrix  R. 
Q  is  formed  as  the  product  of  plane  rotations,  each  which  zero  out  an  entry  in  A.  So 


QA  =  R 


(3.36) 
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The  algorithm  can  be  written  concisely  as 


Kj  =  +  «<+i 

*  =  S,  =  ^ 

0*J  6*,J 

Rows  i  and  i  +  1  are  the  result  of  a  plane  rotation  and  are  given  by 

&ij>  =  +  ^ai+iiP 


(3.37) 

(3.38) 


Ot+ij  —  0  (3.39) 

at+l,P  =:  S*ai,P  d" 

for  p  ^  j.  The  above  is  used  iteratively  until  A  is  reduced  to  triangular  form.  From 
the  above  algorithm  a  data  dependence  graph  can  be  constructed.  This  is  shown  in 
Figure  3.2.  The  arcs  indicate  which  way  data  must  flow  in  order  to  have  the  data 
necessary  to  calculate  the  entry  in  each  node.  This  implies  the  QR  triangularization 
could  be  performed  on  a  two  dimensional  mesh  of  processors.  In  this  case  the  entries 
in  each  node  would  be  calculated  by  a  processor.  This  is  the  implementation  given 
by  Bojanczyk.  Notice  in  Figure  3.2  that  if  a  column  of  the  nodes  is  assigned  to  a 
processor,  then  the  algorithm  can  be  performed  by  a  linear  array  of  processors,  where 
Si  and  Ci  are  the  data  which  must  be  passed  between  processors.  Since  for  the  linear 
array  implementation  both  the  matrix  multiplication  and  the  QR  triangularization 
operate  on  columns  of  the  matrix,  they  can  be  joined  together  very  efficiently.  If  the 
number  of  columns  each  operate  on  are  chosen  the  same,  then  all  the  data  needed  for 
the  QR  factorization  is  resident  in  the  node  after  the  multiplication.  No  communica¬ 
tion  needs  to  take  place.  Next  the  detailed  implementation  for  an  Intel  Hypercube  is 
given. 

Intel  Hypercube  Implementation 

A  set  of  processors  connected  in  a  hypercube  configuration  can  be  formed 
into  a  linear  array.  The  particular  parallel  computer  of  interest  is  an  Intel  Hypercube 
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with  eight  processors.  The  Intel  Hypercube  is  a  multiple  instruction  multiple  data 
computer.  In  this  architecture  each  processor  has  its  own  memory  and  is  connected 
to  its  nearest  neighbor  by  a  communication  line.  The  example  problems  to  be  con¬ 
sidered  consist  of  two  states  and  one  control.  Eight  Walsh  functions  will  be  used  to 
approximate  the  state  and  control.  This  gives  an  M  matrix  which  is  24  x  24  and  an 
augmented  matrix  [K\D]  which  is  24  x  25  on  which  to  perform  a  QR  triangularization. 
Now  solution  on  the  eight  node  Hypercube  is  as  follows: 


•  Form  the  matrix  M  by 


M  = 


INn-{PT®l n)Aa  ~(PT®In)  Af, 
0  0 


•  M  is  sent  to  each  processor  and  each  processor  calculates  three  different  columns 
of  MTM  and  adds  the  appropriate  piece  of  L. 

•  The  d  vector  is  calculated  by 

Di  =  -(/wn  -  ( PT  <8>  In)Aa)T(yec(G')  +  [PT  <g>  In]vec(C )) 

£ 

D2  =  ~-^[PT  ®  ln]\e)T(vec{G.)  +  [PT  <g>  In)vec(C)) 

which  are  matrix  vector  multiplies.  Matrix  vector  multiplies  require  only  T'(n2) 
and  therefore  are  performed  serially. 

•  D  is  augmented  to  K  to  form  [K\D\. 

•  QR  triangularization  and  backsubstitution  are  performed  on  [ifl-D]  yielding  the 
Walsh  coefficients. 


Figure  3.3  summarizes  the  implementation 

The  given  parallel  implementation  is  probably  not  the  fastest  or  the  most 
efficient  possible  implementation.  Because  of  the  nature  of  triangularization,  it  is 
difficult  to  place  the  computational  load  evenly  across  the  linear  array.  The  proces¬ 
sors  which  compute  less  zero  elements  do  the  most  work.  The  implementation  does 
however  show  how  the  method  can  be  matched  to  a  particular  parallel  architecture 


Walsh  Coeffs 


Figure  3.3:  Parallel  Implementation  on  a  Linear  Array 
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or  number  of  processors.  It  also  shows  parallel  implementation  can  be  very  straight¬ 
forward  when  the  problem  has  been  broken  down  into  simple  matrix  operations. 

Although  a  parallel  implementation  for  a  linear  array  of  eight  processors 
has  been  given,  the  potential  parallelizability  of  the  problem  should  be  pointed  out. 
Both  matrix  multiplication  and  QR  triangularization  may  be  implemented  on  a  two 
dimensional  mesh  of  processors.  For  the  problems  outlined  above,  up  to  576  processors 
may  be  used  for  the  matrix  multiplication  and  300  for  the  QR  triangularization.  A 
larger  parallel  computer  which  has  comparable  communication  vs  computation  time 
would  offer  significant  speed-up  over  the  linear  array.  The  next  section  gives  results 
for  a  numerical  example. 


Computational  Results 

The  following  linear  optimal  control  problem  was  implemented  on  the  eight 
node  Intel  Hypercube, 

PROBLEM  3.1: 

Zj  =  x2  xi(0)  =  —4 

x2  =  2xi  —  x2  -f  u  x2(0)  =  4 

J—  [  (2x* -(- xij  +  0.5 u2)dt  (3.41 ) 

Jo 

A  Walsh  series  using  eight  Walsh  functions  is  used  to  approximate  the  state  and 
control.  Since  the  problem  is  linear,  the  cost  functional  of  (3.14)  is  quadratic  and  can 
be  minimized  in  one  iteration.  The  error  in  the  dynamic  constraint  penalty  term  can 
be  made  small  by  setting  epsilon  to  0.0001.  The  plots  in  Figure  3.4  show  the  Walsh 
function  approximations  against  the  continuous  solutions.  Solution  time  on  one  node 
was  2.5  seconds  while  on  all  eight  solution  time  was  0.7  seconds  for  a  speed-up  of  3.5. 

The  linear  optimal  control  problem  has  been  solved  in  many  ways.  It  is  not 
known  for  being  computationally  intensive.  Solutions  to  nonlinear  optimal  control 
problems  on  the  other  hand  axe  considerably  more  complex  and  computationally 
intensive.  A  goal  of  this  research  is  to  solve  such  problems  in  a  highly  parallel  fashion, 
paving  the  way  for  greatly  reduced  solution  times  on  massively  parallel  computers  or 


(3.40) 
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parallel  processing  arrays.  The  parallel  solution  of  the  nonlinear  optimal  control 
problem  is  the  subject  of  the  next  chapter. 


CHAPTER  4 


PARALLEL  SOLUTION  OF  NONLINEAR  OPTIMAL  CONTROL 

PROBLEMS 


This  chapter  shows  how  the  solution  to  the  linear  time- varying  problem  can 
be  used  to  solve  nonlinear  problems.  Two  specific  types  of  nonlinear  problems  are 
addressed:  nonlinear  dynamic  equations  with  quadratic  cost  and  nonlinear  dynamic 
equations  with  quadratic  cost  and  terminal  constraints  on  the  states. 

Optimal  Control  Problems  with  Quadratic  Cost 

The  method  used  to  solve  the  linear  time  varying  optimal  control  problem 
can  be  extended  to  nonlinear  problems.  Quasilinearization  is  one  approach  which  can 
be  used  [17].  For  example,  consider  the  following  nonlinear  problem 

x(t)  =  x(0)  =  x,  (4.1) 

F(x(t),u(f))  =  [  ' {<  x(t),Qx(t)  >  +  <  u(t),Ru(t)  >}dt  (4.2) 

Jo 

where  x(t)  G  /?",«(*)  G  /T".  Let  x*(f)  and  u*(f)  denote  the  optimal  state  and 
control.  Now  choose  some  specific  x(t)  and  u(t)  which  are  close  to  the  optimal  state 
and  control.  These  nominal  values  of  the  state  and  control  will  be  called  x0(t)  and 
u0(t).  Since  x*(t)  and  it *(<)  and  x0 (t)  and  Uo (t)  are  close  together,  then  the  difference 
between  them  can  be  modeled  by  local  linearization.  For  example  let 


6x(t)  =  x*(t)  —  x0(f) 


(4.3) 


Now  loc?l  linearization  is  achieved  by  performing  a  Taylor  series  expansion  of  (4.1) 
and  keeping  only  first  order  terms. 


6x{t)  =  Vxf(x0(t),  u0(t),  t )  8x(t)  +  Vu/(*o(0»  u°(*)> 0  MO  (4-4) 

So  8x(i )  and  6u(t)  can  be  viewed  as  some  perturbation  about  a  nominal.  Now  use 
(4.3)  with  (4.4)  and  get 

(x(t)  -  ±o(0)  =  VJ(i 0(t),  u0(f),  t){x(t)  -  x0 (t))  5 

+ V„/(x0(0, «o(0, 0M0  “  uo(0) 

or  multiplying  out 

*( 0=  Vxf(x0(t),u0(t),t)  x(t)  +  Vuf(x0(t),u0{t),t)  u(t)  ^ 

+(-io(0  ~  Zo(0  -  «o(0) 


which  is  of  the  form 


where 


x(t)  =  A(t)x(t )  +  2?(t)u(t)  +  C(t) 


A(t)  =  V,/(z0(i),t*>(i),i) 

B{t)=  Vu/(zo(0.«o(0>0 
c{t)  =  -/(*o(0»“o(0>0  -  vx/(xo(0»“o(0>0  MO 
-Vu/(xo(<),«o(0>0  uo(0 

This  is  the  form  of  the  linear  time- varying  optin’’!  control  problem  solved  in  Chapter 
3.  Notice  that  if  the  local  linearization  accurately  modeled  the  difference  between  the 
optimal  and  nominal,  if  the  system  was  linear,  then  the  solution  could  be  found  in 
one  step.  For  nonlinear  systems,  however,  the  higher  order  terms  are  ignored.  As  a 
result,  8x(t)  does  not  give  the  change  to  the  optimal,  but  rather  the  change  to  move 
closer  to  the  optimal.  So  once  the  linearized  problem  is  solved  it  is  possible  to  move 
8x(t)  closer  to  the  optimal  solution  by  letting 

x$+1(0  =  xk(t)  k=  1,2,3, .  (4.9) 


As  the  optimal  solution  is  approached,  then  obviously  £x(f)  — *  0.  In  this  way  it  is 
possible  to  solve  the  nonlinear  optimal  control  problem  by  solving  a  series  of  linear 
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time-varying  optimal  control  problems.  The  explicit  quasilinearization  approach  is 
given  as  follows.  Minimize  the  following 

J(e,xk(t),nk(t))  =  /0*'{£  ||  ek(t,e)  ||*  +|  <  xk(t),Qxk(t)  > 

-f|  <  uk(t),Ruk(t)  >}dt 

for  each  k  =  1,2,3, .  where 

ek(t-  c)  =  xk(t)  -  xk.  -  f  Ak{T)x\r)dT  -  /‘  Bk{r)uk{r)dT  -  /“  £*(7-)^  (4.11) 

Jo  JO  JO 

Since  the  problem  is  linearized,  A(t),  B(t),  and  C(t)  are  given  in  equation  (4.7).  After 
each  k  let 

x5+,(()  =  x“(t)  (4,12) 

At  each  iteration,  8x{t)  and  8u(t)  become  smaller  indicating  that  the  series  of  linear 
solutions  is  converging  to  the  nonlinear  solution.  For  computational  purposes  the 
problem  has  converged  when  ||  8x(t)  ||<  e  and  ||  8u[t)  ||<  e  where  epsilon  is  some 
small  number. 

So  the  nonlinear  optimal  control  problem  can  be  solved  by  solution  of  a 
series  of  linear  time  varying  optimal  control  problems,  each  of  which  can  be  solved  in 
parallel.  Computational  results  are  given  in  the  next  section. 


Computational  Results 


The  Raleigh  servo  problem  was  chosen  to  test  the  method.  This  problem  is 
known  to  cause  numerical  instabilities  when  solving  the  dynamical  equations.  Thus 
the  Raleigh  problem  is  very  difficult  to  solve  using  standard  methods.  An  example  of 
the  problem  is  given  as  follows: 

PROBLEM  4.1: 


x‘i  =  Xj  Xi(0)  =  —5 

x-2  =  —  ii  +  1.4  x2  -  0.14  x,  -I-  4  u  xj(0)  =  —5 


V  = 


rO.S 

I  (xj  +  u3)dt 

JO 


(4.13) 

(4.14) 
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Time  (secs) 

Figure  4.1:  Raleigh  Problem  Solution  with  Eight  Walsh  Functions 

To  solve  the  problem,  eight  Walsh  functions  were  used  resulting  in  a  parallel  im¬ 
plementation  as  given  in  Chapter  3.  The  initial  nominal  trajectories  were  set  to  1. 
Epsilon  was  set  to  0.0001.  After  5  iterations,  j|  6x(t)  ||moi<  0.0055  and  the  algorithm 
was  considered  to  have  converged.  The  cost  was  within  .3  %  of  the  optimal  cost 
of  V*  =  17.00.  Figure  4.1  shows  the  control  for  the  first  five  iterations.  Table  4.1 
shows  the  error  function  e(t;e)  and  the  cost  function  J  for  each  iteration.  Figure 
4.2  shows  the  Walsh  approximation  of  the  optimal  control  plotted  against  a  gradient 
based  solution  to  the  epsilon  problem. 

As  in  the  linear  case,  solution  times  are  of  interest.  The  solution  time  on  one 
node  was  13.6  seconds.  On  all  eight  nodes  solution  time  was  about  3.9  seconds,  for  a 
speed-up  of  3.5.  This  speed-up  is  expected  since  the  solution  to  the  nonlinear  problem 
is  just  a  series  of  linear  problems.  The  next  section  examines  how  to  incorporate 
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terminal  constraint  on  the  state  into  the  given  parallel  solution. 


Optimal  Control  Problems  with  Terminal  Constraints 


The  epsilon  method  can  be  used  to  solve  a  variety  of  optimal  control  prob¬ 
lems  [18,  14].  The  last  section  addressed  one  specific  type  of  nonlinear  problem.  The 
Raleigh  problem  has  nonlinear  state  dynamics  and  quadratic  cost.  In  this  type  of 
problem  the  state  is  free  to  take  on  any  value  at  the  final  time.  Some  problems  may 
demand  that  the  state  take  on  a  specific  value  at  the  final  time.  This  section  shows 
how  the  epsilon  method  may  be  used  to  solve  other  forms  of  optimal  control  problems, 
specifically,  those  with  terminal  constraints  on  the  state.  As  before,  the  solution  has 
a  highly  parallel  implementation.  Finally,  to  show  the  effectiveness  of  the  method  the 
well  known  van  der  Pol  problem  is  solved. 

As  shown  in  Frick  [14],  the  optimal  control  problem  can  be  forced  to  meet 
some  final  state  constraint  by  adding  another  penalty  term  which  forces  the  final  state 
to  some  desired  value: 

J(e,x(t),u(t))  =  YcJ''  II  «(*:«)  II’  *+  II  M  II’  +V(*W,u(i))  (4.15) 

As  before, 

e{t\e)  =  x{t)  -  xt  -  /  /(x(T),u(r),r)dr  (4.16) 

Jo 

and 

V(x(t),u(t))  =  /  {<  x(t),Qx(t)  >  +  <  u(t),Ru(t)  >}dt  (4.17) 

Jo 

Approximation  of  the  composite  cost  functional  is  simpler  if  the  new  penalty  term  is 
brought  inside  the  integral.  Now  as  a  quadratic  term  it  is  simplified  similar  to  the 
dynamical  constraint  penalty  term.  The  new  penalty  term  must  force  the  state  to 
some  desired  final  condition  so  let 

p{e)  =  x{tf)  -  Xf  (4.18) 

Here  x/  is  the  desired  final  state.  Since 

*(</)  =  *.  +  f  1  f(x(t),u(t),t)dt 
Jo 


(4.19) 
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then 

p(e)  =  xt+  f(x(t),  u(t),  t)dt  -  x}  (4.20) 

Jo 

Now  consider  the  case  of  linear  time  varying  state  dynamics.  Then 

p(e)  =  x,  +  [  *  A(t)x(t)dt  +  f  1  B[t)u{t)dt  -f  /  C(t)dt  —  Xf  (4.21 ) 

Jo  Jo  Jo 

Note  that  this  is  a  constant  vector  that  can  be  approximated  by 

p{t)  =  m)  (4-22) 


Since  the  vector  to  be  found  is  constant  the  Walsh  approximation  will  be  exact  and 
only  the  coefficients  of  the  zeroth  Walsh  function  will  be  non-zero.  The  cost  function 
J  must  be  put  in  vec  form  so  vec(()  must  be  found. 

Consider  first  the  scalar  case.  Then 

/*'  A(t)x(t)dt  «  f  aV(mT(t)XTdt  (4.23) 

Jo  Jo 


Using  the  orthogonal  property  of  Walsh  functions  this  can  be  written 

f1  aV(t)<ZT(t)XTdt  =  aXT 
Jo 

For  the  case  where  n  >  1 

So  an(t)xi(t)dt  +  •  •  •  +  fo  aln(t)xn(t)dt 
Jo  ani(t)xi(i)dt  + - h  /o  ^nn  (t)xn(t)di 


f  A(t)x(i)dt  = 
Jo 


The  above  matrix  can  be  written 


“ll  •  •  *  <*ln 

’  X?  ' 

®nl  ‘  *  '  ®nn 

(4.24) 


(4.25) 


(4.26) 


Here  the  X  vector  is  not  in  standard  vec  form,  but  it  can  be  rearranged  into  standard 
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vec  form 


(4.27) 


This  matrix  vector  multiply  gives  the  coefficients  of  the  zeroth  Walsh  function,  the 
rest  are  zero.  To  get  a  matrix  of  coefficients  for  the  entire  Walsh  vector,  a  matrix  Ya 
is  formed.  The  first  n  rows  are  obtained  from  (4.26)  and  the  rest  of  the  matrix  is 
zero.  Using  ra  it  is  possible  to  write 


f1  a*(t)*T(t)Xdt  =  Yavec{X) 
Jo 


(4.28) 


The  same  approach  can  be  used  to  approximate  the  integral  with  B(t)u(t).  For  the 
last  integral  a  slightly  different  approach  is  used.  Let 


Then 


<7(t)  =  (7tf(0 


/ 1  C9(t)dt  -  C  f  9(t)dt 
Jo  Jo 


(4.29) 


(4.30) 


The  above  integral  is  a  vector  with  the  first  element  being  a  one  and  the  remaining 
N  are  zero.  The  integral  can  be  written  in  vec  form 


f  CV(t)dt  =  Tcvec(C) 
Jo 


(4.31) 
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Now  vec(()  can  be  written 

vec(()  =  vec(Gt)  +  ravec(X)  +  Tpvec(U)  +  Tcvec(C )  —  vec(Gf)  (4-32) 

To  get  the  system  of  equations  to  be  solved,  the  cost  function  is  minimized  as  before 
with  the  additional  penalty  term  included.  The  result  is  to  add  another  term  to  the 
K  and  D  matrix. 

Kn  =  '-(In.  -  (PT  ®  I.)A.)t(In.  -  (PT  0  /„)A„)  +  (/«  0  (?)  +  ir£ro 
Ku  =  -'-(In.  -  (PT  ®  I.)A.)T([PT  ®  I.]Ae)  + 

*2.  =  -i([pT  ®  i.)Ai>)t(In,  -  (PT  ®  I.)  A„)  +  irjr„  14  33) 

K22  =  \(\PT  ®  /„]Afl)T([/>T  ®  7„]A«)  +  (In  ®  7?)  + 


O.  =  '-(In.  -  (PT  ®  I.)Kf(vec(G.)  +  [PT  ®  /„]v«(G)) 
+;rI{C/  -  G.  -  rcuec(G)} 

Dj  =  -}([PT  ®  /„]A5)T(vec(G.)  +  [/>T  ®  7„]vec(C)) 
+;rJ{G,-G.-i>ec(C)} 

As  done  previously,  the  above  can  be  written  more  concisely  as 


(4.34) 


M  = 


INn-(PT®In)Aa  ~(PT®I n)A0 
0  0 


L  = 


In  ®  Q  0 

0  Iff  ®i  R 


F  = 


rQ  r* 
0  0 


(4.35) 

(4.36) 

(4.37) 


Then 


K  =  -MtM  +  L+-FtF 
e  e 


(4,38) 
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The  extra  term  only  slightly  affects  the  computational  load  of  the  problem  since  F  is 
made  up  of  predominantly  zero  elements.  Since  the  additional  operations  are  simple 
matrix  types,  the  parallelism  of  the  solution  is  not  effected. 

Computational  Results 

Now  a  problem  with  van  der  Pol  oscillator  dynamics  is  solved  on  the  eight 
processor  Hypercube.  The  problem  is  similar  to  the  Raleigh  problem  in  that  it  has  a 
reputation  for  being  difficult  to  solve  as  a  result  of  numerical  instabilities  encountered 
when  solving  the  dynamical  equations. 

PROBLEM  4.2: 

xx  =  x2  *i(0)  =  1  *i(5)  =  --97  ^  39^ 

x2  =  —x\  +  x2  —  x\x2  +  u  x2(0)  =  0  x2(5)  =  —.96 

V  =  J  (x*  4-  x\  +  u2)dt  (4.40) 

Since  the  gamma  matrices  are  nearly  all  zero,  they  can  be  generated  quickly  and  it  is 
not  really  necessary  to  do  a  full  matrix-matrix  multiplication.  The  additional  term 
in  generating  the  D  vector  can  also  be  found  quickly  as  a  result  of  the  sparse  gamma 
matrices.  The  parallel  implementation  is  basically  the  same  as  given  previously, 
except  F  is  added  when  creating  K  after  the  parallel  matrix  multiply  and  an  extra 
term  is  added  to  D  before  it  is  used  in  solution  of  the  system  of  equations. 

To  be  able  to  compare  results  to  previous  solutions,  eight  Walsh  functions 
were  used  to  solve  the  problem.  Epsilon  was  set  to  0.0001  and  the  nominal  trajectories 
were  set  to  unity.  As  in  the  Raleigh  problem,  after  about  five  iterations  the  method 
had  converged.  Table  4.2  gives  statistics  for  the  first  five  iterations.  Figure  4.3  shows 
the  resulting  approximated  optimal  state  and  control.  As  expected,  solution  time 
was  also  similar  to  the  Raleigh  problem.  Solution  time  on  one  node  averaged  14.2 
seconds,  on  eight  nodes  it  averaged  4.1  seconds,  for  a  speed-up  of  around  3.4. 

The  next  chapter  examines  other  orthogonal  functions  which  may  be  used 
effectively  to  solve  the  epsilon  problem.  The  motivation  is  to  find  a  set  of  basis 
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Iteration 

II  HO  llmar 

"MS3DT 

£ 

£ 

J(e,x(t),u(t)) 

1 

1.7326 

0.0307 

0.0034 

12.5862 

2 

0.9912 

0.0025 

0.0001 

3.4257 

3 

0.2668 

0.0040 

0.0002 

4.2447 

4 

0.0057 

0.0041 

0.0002 

4.2490 

5 

0.0004 

0.0041 

0.0002 

4.2490 

Table  4.2:  Iterations  1-5  for  the  van  der  Pol  Solution  with  Eight  Walsh  Functions 


functions  which  might  give  closer,  smoother  approximations  to  the  optimal  while 
using  less  functions.  This  would  mean  a  closer  approximation  with  less  computation. 


CHAPTER  5 


PARALLEL  SOLUTION  OF  NONLINEAR  OPTIMAL  CONTROL 
PROBLEMS  USING  OTHER  BASIS  FUNCTIONS 


Walsh  functions  have  been  applied  to  many  areas  of  science  and  engineering. 
Of  more  specific  interest  here,  they  have  been  applied  to  problems  of  identification, 
analysis,  and  control.  More  recently,  orthogonal  polynomials  have  also  been  employed 
in  optimal  control  methods  [19,  20,  21,  22].  Orthogonal  polynomials  have  the  same  in¬ 
tegration  properties  as  Walsh  functions  and  may  be  used  effectively  as  basis  functions 
to  solve  the  epsilon  problem  when  using  the  Ritz  method. 

In  general,  orthogonal  polynomials  have  the  following  property, 


j  w{t)4>i{t)(f>j{t)dt 
Ja 


0,  j 

ll&ll.  i  =  j 


(5.1) 


A  polynomial  is  orthogonal  with  respect  to  some  weight  function  w(t).  The  weight 
function  for  Legendre  polynomials  is  unity,  similar  to  Walsh  functions.  In  the  follow¬ 
ing  section  some  specific  properties  of  shifted  Legendre  polynomials  are  reviewed. 


Legendre  Polynomials 

Legendre  polynomials  are  defined  on  the  interval  [—1,1].  For  practical  ap¬ 
plication  the  polynomials  are  shifted  to  some  interval  of  interest,  t  6  [0,  </].  The 
shifted  Legendre  polynomials  are  given  by  the  recurrence  relationship  [20] 

(i  +  l)3i+i(t)  =  (2 i  +  \)(~  -  l)si(t)  -  isi.i(t)  i  =  1,2,3,.... 

*/ 


(5.2) 
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1.00 

0.50 


0.00 


-0.50 

-1.00 

0.00  0.20  0.40  0.60  0.80  1.00 

t 

Figure  5.1:  Shifted  Legendre  Polynomials 

with 

s0(t)  =  1 

Figure  5.1  shows  t^e  first  four  shifted  Legendre  polynomials  on  [0,1].  They  are 
orthogonal  on  [0,  t /]  with 

j C  >*  -  ( °;‘  .  (5.3) 

A  function  f(t)  that  is  square  integrable  on  the  interval  [0,  i /]  may  be  expanded  into 
a  shifted  Legendre  series 

/(0- £**(0  (5-4) 

i=0 

Hwang  and  Chen  [20]  showed  that  shifted  Legendre  polynomials  have  inte- 
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gration  properties  similar  to  Walsh  functions,  for  example 

/‘  SN{T)dT  =  HNSN{t)  (5.5) 

J  o 

where  SN{t)  =  [30(t)3i(t)s2(t)  •  •  •  s/sr_i(t)]T.  Thus  Ss(t)  is  a  vector  of  shifted  Legendre 
polynomials  and 

1  1  0  0  •••  0  0  0 

-§  0  |  0  •  ••  0  0  0 

0  — I  Of---  0  0  0 

:  ;  (5.6) 

0  0  0  0  2JV— 3  0  2N-3 

0  0  0  0  •••  0  5^  0 

The  product  of  two  shifted  Legendre  polynomials  is  a  polynomial  and  can 
therefore  be  approximated  by  a  truncated  shifted  Legendre  series.  Hwang  and  Chen 
give  an  approximation  to  the  cross-product  of  two  shifted  Legendre  vectors  as 

50(t)s0(t)  3o(05i(*)  •••  s0(t)sN-i(t) 

si(t)30(t)  ai(t)si(<)  ••• 

SN-i(t)s0(t)  SN-l(t)Si(t)  ■■■  SN~i(t)sN_i(t) 

where 

so(t)s0(t)  =  Zt’o1  2J^9oo  k^t) 

SO(0al(0  =  Zk=0  ^ 901k3k(t ) 

*(0-°(0  =  Efc1 

-»(0-i  (0  =  Efco  2J^9nkSk(t) 

30(t)ss-i{t)  =  'EkJo 

*N-i(i)si(t)  =  T,k=o  2]jfL9N-i,o,k3k(t) 

*N-l(t)sN-l(t)  =  T,k= 0  ^J^9N-l,N-l,kSk{t) 

and 

9i,k  =  [  ’  3t(t)3j(t)Sk(t)dt 

Jo 


(5.8) 


(5.9) 


5(t)5T(t)  = 
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Now  the  same  properties  which  were  defined  for  the  Walsh  functions  are 
defined  for  the  Legendre  polynomials  paving  the  way  for  their  use  as  basis  functions 
for  the  Rayleigh-Ritz  method. 


Rayleigh-Ritz  Method  using  a  Legendre  Polynomial  Basis 


There  are  two  major  differences  between  the  previously  outlined  properties 
of  the  Legendre  polynomials  as  compared  to  the  Walsh  functions.  The  first  is  the  or¬ 
thogonal  property.  For  Walsh  functions  the  result  of  the  orthogonal  property  integral 
is  the  identity  matrix. 

r 1  (  0,  t  ^  j 

f  {  .  (5-10) 

Jo  [  1,  t  =  j 

For  Legendre  polynomials  the  matrix  is  diagonal,  but  it  is  not  the  identity  matrix: 


/■«/ 

/  3i(t)Sj(t)dt 

Jo 


(  0, 

1  aT+T  >  *  =  j 


or  the  result  can  be  written  as  a  matrix 


Ot 


if  0  0 
0  0 
0  0  % 
0  0  0 
0  0  0 


0  0 
0  0 
0  0 
0 


*  =  0, 1,2,  ...N  —  1 


(5.11) 


(5.12) 


Second,  the  cross-product  of  two  vectors  of  the  functions  is  obtained  differently.  For 
Walsh  functions  it  is  obtained  by  simply  multiplying  the  Walsh  functions  out  (see 
Chapter  3).  The  group  property  of  Walsh  functions  guarantees  a  resultant  matrix  in 
which  all  the  resultant  Walsh  functions  are  in  the  same  group  as  the  original  product 
matrices.  For  Legendre  polynomials,  however,  the  cross  product  of  the  two  function 
vectors  is  a  matrix  of  approximations  of  the  product  of  the  individual  polynomial 
multiplications.  This  means  the  result  of  the  product  of  two  Walsh  functions  can  be 
described  exactly  by  another  Walsh  function  whose  subscript  is  no  higher  than  the  two 
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Walsh  functions  being  multiplied.  The  product  of  two  Legendre  polynomials,  however, 
must  be  approximated  by  another  Legendre  polynomial  which  has  the  higher  order 
terms  truncated  to  obtain  a  polynomial  of  the  same  order  as  those  being  multiplied. 
Now  the  solution  using  Walsh  functions  can  be  modified  to  allow  the  use  of  Legendre 
polynomials  and  only  the  above  mentioned  differences  need  to  be  taken  into  account. 

The  composite  cost  functional  (3.5)  in  Chapter  3  must  be  approximated 
with  Legendre  polynomials. 

J=  f\^-ST(t)ETES(t)  +  ]-ST(t)XTQXS(t)  +  \sT{t)UTRUS{t)}dt  (5.13) 

J  0  Z£  4  4 


The  quadratic  terms  in  (5.13)  can  be  simplyfied 

N- 1  N- 1 


ST(t)MS(t )  =  *(0  m«jai(0 

i=0  j= 0 

(5.14) 

and  therefore 

rh  „  n-in-i  t 

/  ST(t)MS(t)dt  =  Y]  mij  /  3i(t)sj(t)dt 

J°  ,= 0  3=0  Jo 

(5.15) 

Since 

[  ’  3i(t)sj(t)dt  =  Oi 

Jo 

(5.16) 

and  if  0/  is  diagonal,  then 

N-  1AT-1  tf  N-l 

H  m'i  /  Si(t)sj(t)dt  =  ^  muon  =  trace(MOi) 

n  w—n  J®  m—n 

(5.17) 

Now  J  can  be  written  in  vec  notation 

J  =  —vec(E)T(Oi<giIn)vec(E)  +  l-vec(X)T(Oi®Q)vec(X)-\-  ]-vec(U)T(Oi® R)vec(U) 

Z  Z 

(5.18) 

The  cost  function  can  now  be  minimized  by  taking  the  gradient  simultaneously  with 
respect  to  both  vec(X)  and  vec(U)  and  setting  them  equal  to  zero.  By  collecting 
terms,  the  result  can  be  written  in  the  same  system  of  equations  format  as  previously 
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given,  where 


Kn  =  \{INn  -  ( PT  ®  In)Aa)T(Oi  ®  7n)(7*n  -  (PT  ®  /„)Aa)  +  (0,  ®  Q) 
Ku  =  -|(7*n  -  (PT  ®  /n)Aa)T(Oi  ®  /n)([PT  ®  /«]A„) 

P«  =  -K(PT  ®  /JA^)T(0,  ®  7n)(7„n  -  ( PT  ®  /„)Aa) 

P22  =  i([pT  ®  7n]A^)T(Oi  ®  /n)([PT  ®  7n]A*)  +  (O,  ®  P) 


(5.19) 


Pi  =  \{lNn  -  (PT  ®  /n)Aa)T(0,  ®  7B)(vec(G.)  +  [PT  ®  7n]uec(C)) 

P2  =  -i([Pr  ®  /n]A^)T(0,  ®  In)(vec(G.)  +  [Pr  ®  In]vec(C)) 

If 

/*„  -  ( PT  ®  /»)A„  -(PT  ®  7n)A^  ■ 

0  0 

and 

_  o,  ®  Q  0 
0  0,  ®  P 

and  now  an  additional  matrix  is  necessary, 


(5.20) 

(5.21) 

(5.22) 


G=  O,  ®  7n  0 
[0  0 

then  P  is 

P  =  +  I 


(5.23) 


(5.24) 


So  the  linear  time  varying  optimal  control  problem  can  be  solved  using  Legendre 
polynomials  as  basis  functions.  As  a  result  of  the  above  development,  the  nonlinear 
optimal  control  problem  can  be  solved  using  Legendre  polynomials.  Only  minor  mod- 
'ications  must  be  made  to  the  parallel  implementation  used  for  the  Walsh  functions. 
The  matrix  G  is  a  diagonal  matrix,  so  it  just  scales  the  columns  of  MT  and  the  full 
three  matrix  multiply  remains  a  two  matrix  multiply. 
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Consider  now  the  addition  of  terminal  constraints  on  the  state.  The  devel¬ 
opment  of  Chapter  4  can  be  followed  with  modification  for  Legendre  polynomials. 


The  penalty  term 

p(e)  =  xt  -f  f  A(t)x(t)dt  +  f  B(t)u(t)dt  +  f  C{t)dt  —  Xf  (5.25) 

Jo  Jo  Jo 

must  be  approximated  with  Legendre  polynomials  in  order  to  find  vec(£).  As  in 
Chapter  4  consider  the  scalar  case  first. 

f  '  A(t)x(t)dt  =  f  aS(t)ST(t)XTdt  (5.26) 

Jo  Jo 

Which  using  the  orthogonal  property  of  Legendre  polynomials  can  be  written 

/*'  aS{t)ST{t)XTdt  =  cxOiXT  (5.27) 

Jo 

Now  for  the  case  where  n  >  1 

{  /o'  an(t)xi(t)dt  +  •  •  •  +  fo1  aln(t)xn(t)dt 

['  A{t)x{t)dt  =  i  (5.28) 

J  0 

So  a.ni(t)xi(t)dt  +  ••■  +  Jo1  ^nn  (t)xn(t)dt  J 
The  above  matrix  can  be  written 


QuOj  OtinOl  X i 

:  :  •  (5.29) 

an\Ol  •••  OnnOl  X% 

This  is  not  in  standard  vec  form,  but  as  shown  in  Chapter  4,  it  can  be  rearranged 
into  vec  notation  and  Ta  can  be  formed.  This  allows  the  integral  to  be  written, 


f'  aS{t)ST(t)XTdt  =  Tavec(X) 
Jo 


(5.30) 


The  same  approach  can  be  used  to  approximate  the  integral  with  B(t)u(t),  and  the 
approximation  for  the  integral  with  C(t)  is  the  same  as  in  chapter  4.  Now  vec(( )  can 
be  written 


vec(()  —  vec(G,)  +  T  avec(X)  +  Tpvec{U )  +  rcvec(C7)  —  vec(Gf ) 


(5.31) 
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The  cost  function  becomes 

J  =  ±vec(E)T(Oi  0  In)vec(E)  +  £vec(0T(O/  0  In)vec(() 
-\-\vec(X)T \Oi  0  Q)vec(X)  +  i[t>ec(C7)T(0i  0  R)vec{V ) 
Minimization  gives  the  system  of  equations 

#11  =  \{.lNn  ~  ( PT  0  7n)Aa)T(0/  0  7„)(7tfn  ~  {PT  0  7n)Aa)  +  (Oi  0  Q) 
+  \Tl{Ot®In)Ta 

#12  =  -}(/*.  -  ( PT  0  In)Aa)T(Ot  0  /n)([PT  0  7n]A*)  +  }rl(0,  ®  /-)r„ 
#21  =  -KfpT  ®  /JA,)r(o,  0  in)(iNn  -  (PT  0  /„)Aa)  +  jrj(o, 0  /n)ra 


(5.32) 


#22  =  f([pT  ®  A.]A*0T(O,  0  /„)([PT  0  7n]A^)  +  (O,  0  72)  +  *rJ(o,  0  7n)r„ 


-Di  =  ;(7w„  -  ( PT  0  7n)Aa)T(0/  0  In)(vec(Gt)  4-  [PT  0  7n]uec(0)) 
+  jrl(0,  0  7„){G/  -  G.  -  Tcvec(C)} 

D?  =  —\{{PT  0  fn]A/3)T(0j  0  7„)(uec(G,)  +  [PT  0  7n]uec(C)) 

+  irj(0,  0  In){Gf  -G.-  Tcvec{C)} 


If 


M  = 


/n„  -  (PT  0  7n)Aa  ~(PT®In)Ap 
0  0 


L  = 


Oi®Q  0 
0  Oj  0  72 


G  = 


0/0  7n  0 
0  0 


F  = 


rQ 

o  o 


then 


#  =  -AfTGM  +  L  +  -FtG# 


(5.33) 


(5.34) 


(5.35) 

(5.36) 

(5.37) 

(5.38) 


Now  both  the  Raleigh  and  van  der  Pol  problems  can  be  solved  using  Legendre  poly¬ 
nomials. 
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Iteration 

II  Sx(t)  ||max 

U-Mll1 

£ 

J(e,x(t),u(t)) 

1 

4.1064 

0.0009 

18.1504 

2 

0.3242 

0.0008 

16.9802 

3 

0.1082 

0.0008 

16.9553 

4 

0.0224 

0.0008 

16.9541 

5 

0.0043 

0.0008 

16.9540 

Table  5.1:  Raleigh  Problem  Solution  with  Eight  Legendre  Polynomials 


Computational  Results 


In  this  section  the  Raleigh  problem  (Problem  4.1)  and  van  der  Pol  problem 
(Problem  4.2)  are  solved  using  Legendre  polynomials.  For  comparison  purposes, 
eight  Legendre  polynomials  are  used.  As  shown  in  the  previous  section  once  the 
product  matrix  S(t)ST(t )  is  determined,  the  only  other  difference  is  the  addition  of 
the  non-identity  orthogonal  matrix  0/.  Since  it  is  a  diagonal  matrix,  full  matrix- 
matrix  multiplications  are  not  necessary. 

Raleigh  Problem  Solution  with  Legendre  Polynomials 

For  this  problem  the  initial  nominal  trajectories  were  set  to  unity.  Epsilon 
was  set  to  0.0001.  As  with  the  Walsh  functions,  after  five  iterations  the  cost  was  within 
.3  %  of  the  optimal  cost,  and  no  noticeable  improvement  in  the  control  was  observed. 
Figure  5.2  shows  the  control  for  the  first  five  iterations.  Table  5.1  shows  the  error 
function  and  the  cost  function  J  for  each  iteration.  Figure  5.3  shows  the  Legendre 
approximation  of  the  optimal  control  plotted  against  the  Walsh  approximation.  The 
solution  time  on  one  node  was  13.6  seconds.  On  all  eight  nodes  the  solution  time  was 
about  3.9  seconds,  for  a  speed-up  of  3.5. 

Actually,  with  the  Legendre  polynomials  it  may  be  possible  to  get  a  very 
close  approximation  to  the  solution  of  the  Raleigh  problem  with  much  less  computa¬ 
tion.  Figure  5.4  shows  the  optimal  control  approximated  by  four  Legendre  polyno¬ 
mials  plotted  against  the  solution  using  four  Walsh  functions. 
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Iteration 

II  £*(i)  Umax 

e 

Mflli! 

e 

xi(5) 

x2(5) 

1 

2.8098 

0.0270 

0.0031 

-.9551 

-.9492 

12.2379 

2 

2.3719 

0.0028 

0.0001 

-.9387 

-.7050 

3.8899 

3 

0.9453 

0.0040 

0.0002 

-.9723 

-.9561 

4.2096 

4 

0.0196 

0.0040 

0.0002 

-.9707 

-.9662 

4.2169 

5 

0.0004 

0.0040 

0.0002 

-.9707 

-.9658 

4.2166 

Table  5.2:  van  der  Pol  Solution  with  Eight  Legendre  Polynomials 


van  der  Pol  Problem  Solution  with  Legendre  Polynomials 

As  with  the  Walsh  functions,  results  for  the  van  der  Pol  problem  (Problem 
4.2)  are  very  similar  to  results  from  the  Raleigh  problem.  The  nominal  trajectories 
were  set  to  unity.  Epsilon  was  set  to  0.0001.  A  solution  was  obtained  in  4-5  itera¬ 
tions.  Table  5.2  gives  the  statistics  of  interest  for  these  iterations.  Figure  5.5  shows 
the  Legendre  approximation  of  the  control  and  states  plotted  against  the  Walsh  ap¬ 
proximation.  Notice  also  from  Table  5.2  that  the  terminal  constraint  is  almost  met 
exactly.  As  expected,  the  solution  times  for  the  case  of  Legendre  polynomials  are  not 
significantly  different  than  for  Walsh  functions.  Solution  time  on  one  node  was  14.3 
seconds  and  on  all  eight  nodes  was  4.2  for  a  speed-up  of  3.4.  Legendre  polynomials 
are  only  one  type  in  the  class  of  orthogonal  polynomials.  It  is  possible  to  use  other 
orthogonal  polynomials  to  solve  the  epsilon  problem.  In  the  next  section  the  shifted 
Chebyshev  polynomials  are  examined. 

Chebyshev  Polynomials  as  Basis  Functions 


It  is  evident  from  Chapter  3  and  5  that  three  key  properties  of  orthogo¬ 
nal  basis  functions  are  used  in  the  solution  of  the  epsilon  problem  using  the  Ritz 
method.  These  properties  result  in  matrices  which  are  used  in  the  solution  to  the 
epsilon  problem.  The  matrices  which  must  be  found  are  the  orthogonal  matrix,  the 
matrix  of  integration  and  the  product  of  two  vectors  of  the  orthogonal  functions.  For 
Walsh  functions,  these  matrices  were  defined  as  Ow,  P,  and  ^'(t)1lI,T(t).  For  Legendre 
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polynomials,  they  were  defined  as  Oi ,  H,  and  S(t)ST(t).  Once  these  matrices  are 
found  for  a  particular  set  of  basis  functions  it  is  straightforward  to  use  the  functions 
as  basis  functions  to  solve  the  epsilon  problem.  This  is  demonstrated  with  Chebyshev 
polynomials  which  are  in  the  class  of  orthogonal  polynomials. 

The  Chebyshev  polynomials  are  defined  between  interval  -1  and  1.  They 
can  be  shifted  onto  the  interval  [0,1]  and  are  given  as 


T0(t)  =  1 

Ti(t)  =  1  —  2t 

T2(t)  =  8ta  -  8t  +  1 

T3(t)  =  -32 13  +  48 13  -  18t  +  1 


(5.39) 


Ti+1{t)  =  {2-4i)Ti(t)-Ti_1(t) 


Multiplication  of  Chebyshev  polynomials  is  not  difficult.  The  cross  product 
of  two  Chebyshev  vectors  is  given  as 

mm  =  i|Tl+i(l)  +  Tv.,1  (5.40) 

Obviously,  to  keep  the  number  of  polynomials  used  as  a  basis  the  same,  the  terms  in 
the  cross  product  matrix  need  to  be  truncated.  For  example  if  N  =  3 


T(t)TT(t ) 


T0(t )  r,(<)  T2(t) 

3*1(0  |[3’2(0  + 3*0(01  §3*1(0 
3i(0  §3*i(0  §3o(0 


where  the  (3, 1),(3,2),  and  (3,3)  terms  have  been  truncated  so  there  are  no  polyno¬ 
mials  of  higher  order  than  T2(t). 

The  integration  matrix  for  shifted  Chebyshev  polynomials  was  derived  by 
Liu  and  Shih  [22],  It  is  given  as  follows, 


fT{t)dt  =  HcT(t ) 
Jo 


(5.41) 
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where  Hc  is  the  integration  matrix 


1 

2 

1 

2 

0 

0 

0 

0 

0 

0 

1 

8 

0 

1 

8 

0 

0 

0 

0 

0 

1 

6 

1 

4 

0 

1 

12 

0 

0 

0 

0 

1 

16 

0 

1 

8 

0 

1 

16 

0 

0 

0 

(5,42) 

1 

2N(N-2) 

0 

0 

0 

0 

1 

4(^-2) 

0 

1 

4AT 

1 

0 

0 

0 

0 

0 

4(JV-1) 

0 

Finally  the  orthogonality  condition  for  shifted  Chebyshev  polynomials  is 


given  as 


[i  mm 

Jo  t-t 2 


I 


0, 

tt/2, 

*■> 


i^j 
i  =  i  7^  o 

i  =  j  =  0 


(5.43) 


The  weight  function  for  Chebyshev  polynomials  is  not  one.  To  remain  consistent  with 
the  parallel  solution  developed  with  Legendre  polynomials,  the  integral  of  the  cross 
product  of  two  Chebyshev  vectors  is  needed.  This  only  needs  to  be  calculated  once, 
so  for  the  given  example  it  is  found  simply  by  numerical  integration.  For  example,  if 


1.0000 

0.0000 

-0.3333 

0.0000 

0.0000 

0.3333 

0.0000 

-0.2000 

-0.3333 

0.0000 

0.4667 

0.0000 

0.0000 

-0.2000 

0.0000 

0.4857 

Raleigh  Problem  Solution  with  Chebyshev  Polynomials 

Using  the  properties  given  in  the  last  section,  the  Raleigh  problem  (Prob¬ 
lem  4.1)  can  be  solved  using  Chebyshev  polynomials.  For  this  problem,  4  Chebyshev 
polynomials  are  used.  Epsilon  is  set  to  0.0001  and  the  nominal  trajectories  are  set  to 
unity.  Table  5.3  gives  the  solution  using  Chebyshev  polynomials.  After  five  iterations 
the  solution  cost  is  very  close  to  the  optimal  cost  of  17.  Figure  5.6  shows  the  Cheby¬ 
shev  approximated  control  plotted  against  Legendre  and  Walsh  approximations.  Just 
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as  in  the  case  of  Legendre  polynomials,  the  Chebyshev  polynomials  give  a  very  close 
approximation  to  the  optimal  solution. 


CHAPTER  6 


NONLINEAR  MINIMUM  TIME  OPTIMAL  CONTROL  PROBLEMS 


In  the  previous  chapters,  only  polynomial  nonlinearities  in  the  dynamical 
equations  were  considered.  In  this  chapter,  some  interesting  computational  properties 
of  Walsh  functions  are  examined  which  allow  a  broader  class  of  nonlinearities  to  be 
dealt  with.  These  properties  of  the  Walsh  functions  also  allow  the  solution  of  complex 
problems  of  a  more  practical  nature,  e.g.,  nonlinear  minimum  time  optimal  control 
problems  which  contain  more  than  simple  polynomial  nonlinearities  in  the  dynamical 
equations.  Results  from  this  chapter  will  also  allow  the  solution  of  problems  which 
involve  inequality  constraints  on  the  state  or  more  typically  the  control.  Solution  of 
problems  with  inequality  constraints  will  be  dealt  with  in  chapter  7. 

Nonlinear  Operations  on  Walsh  Functions 

In  previous  chapters,  polynomial  nonlinearities  in  the  dynamical  equations 
were  dealt  with  by  using  the  multiplication  properties  of  Walsh  functions.  Other 
nonlinear  operations  such  as  the  inverse,  exponential,  or  sin  and  cos  functions  cannot 
be  dealt  with  in  the  same  fashion.  An  alternate  approach  must  be  used. 

Consider  some  function  /(f)  approximated  by  N  Walsh  functions.  The  ap¬ 
proximation  contains  the  value  of  /(f)  at  N  specific  intervals  and  as  such  can  be 
described  by  a  vector  whose  N  elements  give  the  value  of  /(f)  at  each  of  the  N  inter¬ 
vals.  This  is  demonstrated  by  considering  the  Walsh  function  <^2(f)  (See  Figure  2.2). 
This  Walsh  function  can  be  completely  described  by  the  vector 
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where  each  of  the  elements  gives  the  value  of  fa (t)  on  each  of  its  four  different  intervals. 
As  another  example,  consider  fa(t).  It  is  described  by 

P3  =  [  1  -1  -1  1  ] 

It  is  possible  to  describe  first  4  Walsh  functions  by  a  vector  with  4  elements,  the  first 
8  Walsh  functions  by  vectors  with  8  elements  and  so  on.  In  this  manner,  the  vector 
\?( t )  can  be  described  by  a  matrix, 

fait)  Po 

fa  (t)  _  Pi 

J  Pn 

And,  just  as  it  is  possible  to  approximate  a  function  by 

fit)  »  F9( t)  (6.2) 

the  alternate  form  can  be  used 

fit)  as  FV  (6.3) 

where  V  now  is  a  matrix.  The  following  example  will  help  clarify  the  process. 
EXAMPLE: 

Let  f(t)  =  t  and  N  =  4,  then  /(<)  can  be  approximated  by  a  Walsh  series 

fit)  «  F«( t) 


where 

Thus 


fit)  KFt  =  FV=[  .125  .375  .625  .875 


and  Ft  describes  the  Walsh  series  approximation  for  fit). 

The  operation  just  described  amounts  to  an  efficient  numerical  method  of 
obtaining  the  actual  Walsh  series  from  the  Walsh  coefficients.  It  also  gives  a  technique 
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for  getting  the  Walsh  coefficients  from  the  Walsh  series.  Since  the  Walsh  functions 
are  orthogonal,  then  V  is  invertable  and  it  is  possible  to  write 

FtV~x  =  F  (6.4) 

where  Ft  is  some  Walsh  approximation  to  a  time  function  given  in  the  vector  form  and 
F  then  is  the  matrix  of  Walsh  coefficients  of  the  Walsh  series.  By  use  of  the  V  matrix  it 
is  possible  to  generate  Walsh  functions  from  Walsh  coefficients  and  Walsh  coefficients 
from  Walsh  functions.  Performing  nonlinear  operations  on  Walsh  functions  described 
by  a  set  of  Walsh  coefficients  now  becomes  straightforward.  Nonlinear  operations 
cannot  be  performed  directly  on  a  set  of  Walsh  coefficients,  but  they  can  be  performed 
elementwise  on  a  Walsh  series  in  the  vector  form.  An  example  is  given  to  clarify  this 
point. 

EXAMPLE: 

It  is  desired  to  find  g(t)  =  sin (/(f))  with  only  the  Walsh  coefficients  of  /(<) 
are  known.  In  other  words  , 

m  «  pm 

and  F  is  known.  It  is  desired  to  find  G  where 

g(t)  «  G^(t)  =  sin (F'P(t)) 

This  is  straightforward  using  the  V  matrix.  The  above  is  equivalent  to 

sin(FT’)  =  GV  =  Gt 


so 

G  =  sin  {FV)V~X 

Let  N  =  4  and  assume 

F=  [  1.5  -.75  -.375  0  ] 

which  gives  the  Walsh  approximation  for  f(t)  =  3 1.  From  F,  the  Walsh  coefficients 
for  sin(3t)  are  easily  found 


G  =  sm(FV)V~1  =  [  .679  -.045  -.019  -.249 
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Figure  6.1  shows  the  Walsh  approximation  using  G  against  the  actual  function  sin(3i). 
Other  nonlinear  functions  can  be  performed  on  Walsh  coefficients  in  the  same  manner. 
By  use  of  the  V  matrix,  it  is  possible  to  solve  nonlinear  optimal  control  problems  which 
have  many  types  of  nonlinearities  other  than  just  polynomial  nonlinearities.  In  the 
next  section  this  is  demonstrated  by  numerical  solution  of  a  nonlinear  minimum  time 
optimal  control  problem.  This  problem  has  other  than  polynomial  nonlinearities  in 
the  system  dynamical  equations  and  is  intended  to  demonstrate  the  practical  value 
of  the  method. 


Solution  of  a  Nonlinear  Minimum  Time  Optimal  Control  Problem 

In  the  previous  section,  the  tools  were  developed  to  allow  the  parallel  solu¬ 
tion  of  nonlinear  optimal  control  problems  with  other  than  polynomial  nonlinearities. 
This  will  be  demonstrated  by  solution  of  a  nonlinear  minimum  time  optimal  control 
problem.  The  particular  problem  of  interest  involves  determination  of  the  control 
which  generates  a  minimum  time  transfer  trajectory  from  Earth  orbit  to  Mars  orbit. 
The  following  problem  is  a  variation  of  a  problem  solved  by  Taylor  [23].  The  equa¬ 
tions  of  motion  which  govern  a  spacecraft  with  constant  thrust  throughout  the  time 
of  flight  axe  given  as 

PROBLEM:  6.1: 


Xi 

x2 

x3 


x 2 

x\  n  T  sin  0 

Xj  xj  ma  -f  mt 

x2x3  T  cos  9 

x\  m„  +  mt 


(6.5) 


where  Xj  is  radial  position,  x2  is  radial  velocity,  x3  is  tangential  velocity,  T  is  thrust, 
6  is  thrust  angle,  /z  is  the  gravitational  constant,  m„  is  the  initial  mass,  m  is  mass 
flow  rate,  and  t  is  time.  The  problem  is  normalized  to  astronomical  units.  The  initial 
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values  for  the  state  variables  are 

*i(0)  =  1.0 

*a(0)  =  0.0  (6.6) 

a3(0)  =  1.0 

The  above  initial  conditions  define  an  initial  orbit  which  is  the  same  as  Earth’s  orbit. 
The  final  values  for  the  state  variables  must  be 

xu  =  1.525 

*2(/  =  0.0  (6.7) 

x3,  =s  0.8098 

The  final  conditions  define  a  final  orbit  which  is  the  same  as  Mars’s  orbit.  Parameter 
values  are 

H  =  1.000 

T  =  0.1405 
m0  —  1.000 
m  =  0.0749 

and  1  time  unit  equals  58.18  days.  Finally,  the  cost  function  to  be  minimized  is 

J  =  jf*'  1  dt  (6.8) 

This  defines  the  problem  as  a  minimum  time  problem.  In  other  words,  the  optimal 
time  ty  must  be  found  which  minimizes  J.  The  problem  just  given  can  be  summarized 
as  follows:  given  the  equations  of  motion,  find  the  thrust  angle  trajectory  which 
transfers  a  spacecraft  from  Earth  orbit  to  Mars  orbit  in  the  minimum  time. 

Solution  by  Steepest  Descent 

To  obtain  a  reference  for  later  comparisons,  the  problem  is  first  solved  using 
the  method  of  steepest  descent.  Steepest  descent  is  basically  a  method  for  solving 
the  two  point  boundary  value  problem  resulting  from  (2.15)  thru  (2.19)  in  Chapter  2. 
The  form  of  these  equations  naturally  lend  themselves  to  solution  by  steepest  descent. 
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First  the  TPBVP  which  describes  the  Earth  to  Mars  problem  is  given.  The 
equations  of  motion  are  repeated  for  convenience. 


Xi 

x2 


=  x2 


X3  = 


xZ  u 

---=*  + 
Xi  x\ 

X2X3 


T  sinu 
m0  +  mi 
T  cos  u 


x\  m0  +  mt 


The  Hamiltonian  is 


H  —  if  - f  A1X2  +  A2(— - j  + 


X3  M  ,  rsinu  )  |  W 
Xi  xf  m0  +  ihtJ 


x2x3  T  cos  u  . 
xi  m0  +  mt 


From  this  the  costate  equations  can  be  written 

„2 


Finally 


Ai 

Aj 

A3 

dH 

du 


X  j  X j 

=  — Ai  +  [— ]A3 

Xl 

=  -[— ]A2  +  [^]A3 

Xj  Xi 


[Aj  cos  u  —  A3  sin  u 


m0  +  mt 


(6.9) 


(6.10) 


(6.11) 


(6.12) 


Note  that  the  terminal  conditions  can  be  met  explicitly,  or  the  following  cost  function 
can  be  used 


J  =  Jq  1 1  dt  +  5(ii(</)  -  I,,)’  +  S(xa(tf)  -  x3 t)f)2  +  S(x3(tj)  -  x3t/)2  (6.13) 

which  forces  the  states  to  be  very  close  to  the  required  terminal  state  values  if  5  is 
chosen  very  large.  In  this  case,  the  boundary  conditions  are 


xi(0)  =  1.0 
x2(0)  =  0.0 
x3(0)  =  1.0 


M*/)  =  -xllf) 

A2  (tf)  =  s(x3(if)-x  3t/) 
A 3(t/)  =  s(x3 (if)  -  x3tf) 


(6.14) 


(6.15) 
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Now  a  control  trajectory  u{i )  is  chosen  which  is  close  to  the  optimal  control  trajectory 
u*(t).  The  state  equations  can  be  integrated  using  this  u (t)  and  the  final  conditions 
on  the  state  axe  used  to  find  the  final  conditions  on  the  costate  using  (6.14).  The 
costate  equations  are  then  integrated  backwards.  Now  the  costate  trajectories  can  be 
used  to  determine  If  ^  =  0  then  u(t)  =  u*(t),  otherwise  u(t)  can  be  corrected 

by 


ui+1(y)  =  Ui(t)  -  K  —  (6.16) 

where  K  <  1  and  is  a  small  number  chosen  for  best  convergence.  The  new  u(t) 
replaces  the  old  one  and  the  process  is  repeated  until 


(6.17) 


where  e  is  some  arbitrarily  small  number.  This  is  the  method  of  steepest  descent. 

In  minimum  time  problems,  the  cost  function  must  also  be  minimized  with 
respect  to  the  final  time  tf.  Since  the  cost  function  is  a  function  of  £/,  one  way  to 
minimize  the  cost  function  is  to  solve  a  series  of  steepest  descent  problems.  For  each 
consecutive  problem  choose  the  new  t/  by  some  line  search  optimization  technique, 
for  example  quadratic  curve  fitting  (See  Luenberger  [24]).  This  technique  was  used 
to  solve  the  Earth  to  Mars  minimum  time  problem.  The  starting  control  was  chosen 
to  be  1.2  radians  for  the  first  half  of  flight  and  —1.2  radians  for  the  second  half  of 
flight.  For  S  =  10000,  the  solution  is 


tf  =  3.9515 
j final  =  3.9589 

Figures  6. 2-6. 5  show  the  optimal  control  and  states. 


Solution  by  the  Epsilon  Method 

The  integral  epsilon  formulation  of  the  Earth  to  Mars  minimum  time  prob¬ 
lem  is  given  as 

J(e,x(t),u(t))  =  -{£'  »  II’  dt+  ||  p(c)  II’)  +  J‘‘  Idt  (6.18) 
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as  before 


and 


i{t\e)  =  x(t)  -  xt  -  [  /(®(t),u(t),t)<£7 
Jo 

ftj 

p(e)  =  x,+  f(x(t),u(i),t)dt  —  Xj 
Jo 


(6.19) 


(6.20) 


This  problem  can  be  solved  as  shown  in  Chapters  3  and  4.  First  the  problem  must 
be  linearized 


x(t)  =  j4(t)x(t)  4  B(t)u(t)  +  C(t) 


where 


m  = 


i  o 


+  0  2 


lo  1 
xlox3o 
*lo 


x3c 
XU 

x3o  _  £Jj 
lie  *lo  J 


(6.21) 


(6.22) 


B(l)  = 


cm 


a(t)  cos  ua 
— a(t)  sin  u0 

0 


(6.23) 


4-  a(t )  sin  ua  —  ( a(t )  cos  uo)u0 


a(t)  cos  ua  +  (a(t)  sin  u0)u0 
Here  xlol  X2„,  X30,  and  ua  are  some  nominal  trajectories  and 

T 


(6.24) 


a(t)  = 


m0  +  rht 


The  Walsh  coefhcients  needed  to  approximate  the  above  matrices  can  be  found  using 
the  V  matrix  once  the  initial  nominal  trajectories  have  been  chosen.  Once  the  ma¬ 
trices  have  been  approximated  in  terms  of  the  Walsh  coefficients,  the  problem  can  be 
solved  as  before  in  Chapter  4.  To  find  the  optimal  tj,  some  line  search  optimization 
is  performed  as  done  with  the  steepest  descent  technique.  The  problem  solution  is 
summarized  as  follows: 


Pick  initial  nominal  trajectories. 


•  Pick  an  initial  if. 
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Walsh  Functions 

J final 

4 

4.018 

4.042 

8 

3.949 

3.975 

16 

3.921 

3.946 

Table  6.1:  Minimum  Time  versus  Number  of  Walsh  Functions 

•  Solve  the  linearized  problem  iteratively  until 

||  2/  XQ  ||mijt<  Cl 

•  Change  tj  using  some  optimal  search  scheme  and  repeat  until 

ii  4+i  - 1)  n< 

where  ei  and  62  are  some  small  numbers. 

The  above  scheme  was  used  to  solve  the  Earth  to  Mars  problem  using  4,8 
and  16  Walsh  functions.  The  nominal  trajectories  were  chosen  as  follows: 

xu  =  1.0 
x2o  =  0.0 
x3e  =  1-0 
u0  =  .5 

Table  6.1  gives  the  minimum  time  results  as  a  function  of  the  number  of  Walsh 
functions.  Results  of  a  single  epsilon  problem  using  eight  Walsh  functions,  the  min¬ 
imum  time,  and  the  given  initial  nominal  trajectories  are  summarized  in  table  6.2. 
As  expected,  speed  of  convergence  is  dependent  on  the  initial  pick  of  the  nominal 
trajectories.  Since  tto  is  not  very  close  to  u*(t),  it  takes  a  number  of  iterations  to 
converge.  It  does  show,  however,  that  the  method  converges  even  for  relatively  poor 
initial  guesses.  To  test  the  actual  effectiveness  of  the  method,  the  approximated  con¬ 
trol  was  used  to  drive  the  nonlinear  equations  of  motion  (6.9)  from  the  given  initial 
conditions.  The  results  of  the  simulation  are  plotted  against  the  solution  by  steepest 
descent  in  Figures  6.2  thru  6.5.  Apparently  the  method  gives  very  good  results  even 
for  as  few  as  eight  Walsh  functions. 


Thrust  Angle  u(t) 
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1.40 

1.00 

0.60 


Epsilon  Method 
Steepest  Descent 


0.20 

0.20 


-0.60 

-1.00 


-1.40  1 - 1 - 1 - 1 - 1 - ' - 1 - 1 - 1 

0.00  0.50  1.00  1.50  2.00  2.50  3.00  3.50  4.00 


Time 


Figure  6.2:  Optimal  Thrust  Angle  -  Minimum  Time  Solution 


Radial  Position 
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Figure  6.3:  Radial  Position  -  Minimum  Time  Solution 


Radial  Velocity 
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Time 


Figure  6.4:  Radial  Velocity  -  Minimum  Time  Solution 


Tangential  Velocity 


1.0  5 

1.00 

0.95 

0.90 

0.85 

0.80 

0.75 
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Although  the  minimum  time  example  was  solved  on  a  serial  machine,  im¬ 
plementation  of  the  problem  in  parallel  poses  no  new  problems  and  speed-ups  similar 
to  previous  parallel  examples  can  be  expected. 


CHAPTER  7 


OPTIMAL  CONTROL  PROBLEMS  WITH  INEQUALITY 

CONSTRAINTS 


In  many  optimal  control  problems,  the  state  or  control  must  be  constrained 
to  stay  within  certain  limits.  Consider  the  following  problem 

*(*)  =  /(*(0ft*(0i0  x(t0)  =  x0  (7.1) 

V(*(0.u(0)“  (7.2) 

Jo 

and 

l«(<)|<c  (7.3) 

where  c  is  some  constant.  In  this  case,  it  is  desired  to  find  the  control  which  minimizes 
V  while  meeting  the  dynamic  constraint  of  (7.1)  and  the  inequality  constraint  of 
(7.3).  The  control  is  constrained  to  remain  within  some  specified  limits,  so  this 
is  a  constrained  dynamic  optimization  problem.  Contrary  to  the  previously  solved 
problems,  there  are  two  constraints  which  must  be  taken  into  account:  the  equality 
constraint  of  the  dynamical  equations  and  the  inequality  constraint  on  the  control. 
The  linear  time- varying  optimal  control  problem  with  no  inequality  constraints  was 
solved  in  chapter  3.  By  conversion  into  the  integral  epsilon  problem  and  solution 
using  the  Rayleigh-Ritz  method,  the  dynamic  optimization  problem  was  converted 
into  a  static  optimization  with  the  following  cost  function: 

J  =  -vec(E)Tvec(E)  +  ]-vec(X)T(lN®Q)vec(X)+]-vec(U)T(Iff<g)R)vec(U)  (7.4) 

M  A 


where 


vec(E)  =  { vec(X )  —  vec(G,)  —  \{PT  0  7n)Aa]vec(A’) 
-((PT  8  U)he)vec{U)  -  (PT  8  In)vec(C)} 


(7.5) 
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This  cost  function  is  minimized  simultaneously  with  respect  to  both  the  Walsh  coeffi¬ 
cients  for  the  state  vec(X)  and  the  Walsh  coefficients  for  the  control  vec{U).  If  there 
are  inequality  constraints  on  the  state  or  the  control,  then  the  above  cost  function 
must  be  minimized  while  meeting  the  inequality  constraint.  Since  the  cost  function 
is  minimized  with  respect  to  the  Walsh  coefficients,  the  inequality  constraint  must 
be  formulated  in  terms  of  the  Walsh  coefficients.  Once  this  is  done  the  problem  be¬ 
comes  a  static  constrained  optimization  problem.  Similar  to  the  epsilon  method,  this 
problem  can  be  solved  by  adjoining  the  inequality  constraint  to  the  cost  function  as 
a  penalty  term  and  then  minimizing  with  respect  to  the  Walsh  coefficients.  This  es¬ 
sentially  converts  the  static  constrained  optimization  problem  into  an  unconstrained 
problem.  The  key  then  is  to  be  able  to  formulate  the  inequality  constraint  as  given 
in  (7.3)  into  an  equivalent  inequality  constraint  on  the  Walsh  coefficients.  The  next 
section  shows  how  an  inequality  constraint  on  the  control  or  state  can  be  formulated 
in  terms  of  the  Walsh  coefficients. 


Formulating  Walsh  Coefficient  Inequality  Constraints 


Consider  an  optimal  control  problem  with  the  following  inequality  con¬ 


straint, 


I  “(<)  |<  C 


(7,6) 


where  c  is  some  constant.  This  inequality  constraint  must  be  converted  into  a  con¬ 
straint  on  the  Walsh  coefficients  of  the  control.  The  inequality  constraint  must  be 
given  in  terms  of  the  Walsh  coefficients.  This  can  be  done  using  the  V  matrix  given 
in  chapter  6.  For  example,  c  can  be  described  by  a  vector  of  N  elements  of  value  c. 
Let 


c  =  C,=  [c?  c>  ...  c"-'] 

where  Ct  is  a  vector  of  N  values  for  which  cj  =  c.  In  the  same  manner 


u(t)  =  ut=  [u°  U ]  •••  tx?-1  ‘ 
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Now  the  inequality  constraint  can  be  written 

K  \<c\  i  =  0,1,2,. . .  ,N  —  1 

Since 

Ut  =  UV  =  Up,  i  =  0,1,2,...  ,N  —  1  (7.7) 

the  inequality  can  be  written 

I  UPi  |<c*  i  =  0,1,...  ,N  —  1  (7.8) 

In  (7.8),  pi  are  the  vectors  which  describe  the  Walsh  functions  as  shown  in  (6.1).  The 
next  section  gives  one  method  which  may  be  used  to  solve  the  constrained  optimiza¬ 
tion  problem  given  by  (7.4)  and  (7.8). 


Solution  of  Optimal  Control  Problems  with  Inequality  Constraints 

In  a  manner  similar  to  the  epsilon  problem,  the  constrained  optimization 
problem  given  by  (7.4)  and  (7.8)  can  be  converted  into  an  unconstrained  problem 
by  adjoining  the  inequality  constraints  to  the  cost  function  in  the  form  of  a  penalty 
term.  The  following  penalty  term  may  be  used, 

||  max(0,|  UPi  |  -c*,.)  ||2  (7.9) 

In  this  case,  only  the  constraints  which  are  violated  will  appear  in  the  penalty  term. 
A  good  approach  then  is  to  determine  at  each  iteration  the  constraints  which  are 
active  and  form  a  matrix  V  and  C[  which  consist  of  only  the  active  constraints.  This 
matrix  and  vector  are  then  used  to  form  the  penalty  term  which  must  be  in  vec  form. 

i  II  ®  WIO  -  »««<?,')  II1  (7.10) 


This  term  is  added  to  the  cost  function  of  equation  (7.4).  Let 

vec(S)  =  (P'T  ®  In)vec(U)  —  vec(C,t') 


(7.11) 
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then  the  cost  function  becomes 

J  =  £vec(E)Tvec(E)  +  \vec(X)T(Iff  ®  Q)vec(X ) 

+  ^vec(U)T(Ift  ®  R)vec(U)  -f  j^vec(S)T vec(£) 

Now  the  cost  function  can  be  minimized  as  before  and  the  system  of  equations 


Xu  K\  2 
X 2i  K<n 


vec*(X) 

vec*{U) 


D i 

D2 


is  obtained  by  letting 

M  = 

and 

and 

Then 


INn  -  ( PT  ®  /„) Aq  -(Pt  ®  In)A0 

0  0 


L  = 


In  ®  Q  0 
0  In  ®  R 


F  = 


0  0 
0  V”®In 


K  =  -MtM  +  L+  -FtF 

z  7 


and  a  term  must  also  be  added  to  D 


(7.13) 


(7.14) 


(7.15) 


(7.16) 


(7.17) 


D2  =  D2  +  \vrr  ®  I„)Tvec(C[)  (7.18) 

In  (7.18),  D3  is  the  same  as  developed  in  Chapter  3.  To  incorporate  inequality 
constraints,  at  each  iteration  the  active  constraints  are  determined.  This  involves 
checking  N  intervals  of  the  state  or  control.  If  the  value  of  the  state  or  control 
in  an  interval  exceeds  c\,  that  constraint  must  become  active.  Once  all  the  active 
constraints  are  determined,  the  cost  function  is  minimized  by  solving  the  system 
of  equations  given  by  including  the  modifications  above.  The  next  section  gives  a 
numerical  example  to  demonstrate  the  method. 
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Computational  Results 

To  demonstrate  the  method  given  in  the  previous  section,  following  problem 
was  solved. 

PROBLEM  7.1: 

x\  =  X2  3i(0)  =  0 

35  2  =  —  Xi  +  X2  —  xfx2  +  U  3:2(0)  =  1 

V  =  f\x  \  +  x\  +  u2)dt  (7.20) 

Jo 

This  is  a  van  der  Pol  problem  with  no  terminal  constraints.  Now  Problem  7.2  is 
obtained  by  adding  an  inequality  constraint  to  the  control. 

PROBLEM  7.2: 

35*1  =  X2  2:i(0)  =  0 

X2  =  —xx  +  3:2  —  xjx2  +  u  x2(0)  =  1 

V  =  (\x\  +  x\  +  u2)dt  (7.22) 

JO 

I  u(t)  \<  -8  (7.23) 

To  solve  both  problems,  eight  Walsh  functions  were  used  and  epsilon  and  gamma 
were  set  to  0.0001.  Figure  7.1  shows  the  optimal  control  and  state  for  Problem  7.1. 
Figure  7.2  shows  the  solution  when  the  control  is  constrained  to  stay  between  —.8 
and  .8 

Although  the  above  example  was  solved  on  a  serial  computer,  implementa¬ 
tion  of  the  problem  in  parallel  would  involve  only  minor  modifications  to  the  parallel 
implementation  given  in  Chapter  3.  The  only  addition  would  be  the  determination 
of  the  active  constraints.  Since  each  of  the  N  intervals  are  independent  of  each  other 
this  may  be  done  in  parallel.  As  shown  with  the  addition  of  the  penalty  term  for 
terminal  constraints  in  Chapter  4,  addition  of  another  penalty  term  will  not  effect 
parallelization  of  the  problem. 


(7.21) 


(7.19) 


CHAPTER  8 


CONCLUSIONS 


This  research  has  developed  and  implemented  a  parallel  solution  method 
for  nonlinear  optimal  control  problems.  Balakrishnan’s  epsilon  method  converts  the 
two  point  boundary  value  problem  into  a  nonlinear  programming  problem  thus  avoid¬ 
ing  the  serial  bottlenecks  inherent  in  methods  which  solve  the  dynamical  equations. 
The  solution  to  the  nonlinear  programming  problem  can  be  found  using  simple  ma¬ 
trix  operations  which  have  parallel  implementations  with  the  desirable  properties  of 
flexibility  and  scalability. 

To  test  the  method,  a  variety  of  nonlinear  optimal  control  problems  were 
solved.  The  Raleigh  problem  represents  the  case  of  nonlinear  dynamical  equations 
with  a  quadratic  cost  and  the  van  der  Pol  problem  represents  the  case  of  nonlinear 
dynamical  equations  with  a  quadratic  cost  and  terminal  constraints  on  the  state. 
Using  eight  Walsh  functions  or  Legendre  polynomials  to  approximate  the  state  and 
control,  the  problems  could  be  solved  in  four  to  five  iterations  from  relatively  poor 
initial  guesses  for  the  nominal  trajectories.  On  an  eight  node  Intel  Hypercube,  solution 
time  was  4-5  seconds  and  speed-up  on  the  order  of  3. 4-3.5. 

The  Raleigh  problem  and  van  der  Pol  problem  contain  only  polynomial 
nonlinearities  in  the  system  dynamical  equations.  It  is  simple  to  perform  complex 
nonlineax  operations  on  Walsh  functions,  making  them  ideal  for  the  solution  of  prob¬ 
lems  with  nonlinearities  which  are  more  complex  than  just  polynomial  nonlinearities. 
This  was  demonstrated  by  solution  of  a  nonlineax  minimum  time  optimal  control 
problem  and  nonlinear  optimal  control  problem  with  an  inequality  constraint  on  the 
control.  Results  for  the  nonlinear  minimum  time  problem  were  compared  to  results 
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obtained  by  solution  of  the  same  problem  using  the  technique  of  steepest  descent. 
Results  show  the  method  gives  accurate  results  and  is,  therefore,  a  feasible  technique 
for  the  rapid  solution  of  complex  nonlinear  control  problems. 

Since  it  is  not  apparent  how  to  deal  with  other  than  polynomial  nonlinear¬ 
ities  when  using  orthogonal  polynomials,  it  would  appear  that  the  Walsh  functions 
have  a  broader  range  of  applications.  They  would  also  be  more  effective  when  con¬ 
sidering  problems  with  discontinuities  in  the  state  or  control.  On  the  other,  hand 
preliminary  results  show  that  for  problems  with  continuous  state  and  control  a  good 
approximation  may  be  obtained  using  very  few  orthogonal  polynomials.  As  a  result, 
the  best  basis  functions  to  solve  a  problem  will  be  a  function  of  the  particular  problem 
to  be  solved. 

Finally,  it  should  be  pointed  out  that  the  parallel  implementations  given  in 
this  thesis  were  geared  toward  a  specific  parallel  computer.  Matrix  operations  have 
a  variety  of  parallel  implementations.  As  a  result,  it  is  possible  to  solve  the  optimal 
control  problem  on  a  variety  of  parallel  computers  or  parallel  processing  arrays.  In 
addition  to  this  implementation  flexibility,  potential  parallelization  is  very  high.  For 
example  using  eight  Walsh  functions  up  to  560  processors  might  be  used  to  solve  the 
minimum  time  problem  given  in  Chapter  6.  The  method  offers  the  potential  for  very 
fast  solution  times  of  complex  nonlinear  optimal  control  problems. 

Future  Research  Areas 

This  research  paves  the  way  to  many  additional  research  areas.  Only  a  small 
parallel  machine  was  available  to  solve  the  example  problems.  It  would  be  a  powerful 
demonstration  to  solve  a  complex  nonlinear  optimal  control  problem,  such  as  the 
minimum  time  problem,  on  a  very  large  parallel  computer  (p  >  500).  The  goal  may 
be  solution  times  approaching  real-time. 

Inherent  in  the  given  parallel  solution  technique  is  a  highly  parallel  solu¬ 
tion  technique  for  nonlinear  ordinary  differential  equations.  There  are  other  areas 
of  estimation  and  control  in  which  a  parallel  solution  method  for  solving  ordinary 
differential  equations  would  be  useful.  For  example,  it  may  allow  a  highly  parallel 
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solution  method  for  extended  Kalman  filters  or  parallel  solution  of  Riccati  equations. 

Frick  showed  that  Walsh  functions  could  be  used  to  approximate  noise  pro¬ 
cesses  [25].  It  would  be  interesting  to  investigate  parallel  solution  methods  for  stochas¬ 
tic  optimal  control  problems  using  these  ideas. 

Finally,  the  idea  of  approximating  functions  with  a  set  of  basis  functions 
to  formulate  a  parallel  solution  method  is  a  useful  one.  It  is  possible  the  method  of 
approximating  functions  may  be  used  to  reformulate  other  problems  which  seem  very 
serial  in  nature  into  problems  which  are  highly  parallel  in  nature. 
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Appendix  A 
Notation 


This  appendix  reviews  some  very  powerful  systems  notation  as  given  from 
a  paper  by  Brewer  [15].  The  notation  is  used  in  the  solution  of  the  epsilon  problem 
using  the  Ritz  method  with  some  orthogonal  functions  as  a  basis. 

Let  A  be  a  matrix  of  dimension  p  x  q,  B  is  dimension  sxt  and  D  is  dimension 
q  x  s.  Then  vec  of  A  is  a  vector  of  dimension  pq  x  1  defined  as 


(A-l) 


Here  the  A.\  notation  refers  to  the  first  column  of  A  and  A.2  refers  to  the  second  and 
so  on.  So  the  uec(A)  amounts  to  stacking  the  columns  of  A  on  top  of  each  other  into 
a  large  vector. 

The  symbol  ®  denotes  the  Kronecker  product.  It  is  defined  as: 


So  A  ®  B  is  a  ps  x  qt  matrix. 

Finally,  it  can  be  shown  that 


vec(ADB)  =  ( Bt  ®  A)vec(D ) 


(A3) 
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Using  this  identity  the  vec  of  the  product  of  two  matrices  can  be  found. 

vec(AD)  =  vec(ADIt )  =  ( I ,  ®  A)vec(D) 
vec(AD)  =  vec(IpAD)  =  ( DT  ®  Ip)vec(A) 
vec(AD)  =  vec(AIqD)  =  (DT  ®  A)vec(Iq ) 


(A4) 


